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SUMMARY 

During the past 6 months our effort has been concentrated on the following 
3 general areas: 

(1) Continued development of realistic, finite element modeling of plate 
interactions and associated deformation in the Eastern Mediterranean, 

(2) Neotectonic field investigations of seismic faulting along the active fault 
systems in Turkey with emphasis on identifying seismic gaps along the 
North Anatolian faulty 

(3) Establishment of a CPS regional monitoring network in the zone of ongoing 
continental collision in eastern Turkey (supported in part by NSF). 

A short summary is given below. More detailed descriptions of each of these 
aspects of our research are given in the attached Appendices. 

Finite Element Modeling (Appendix l) 

The purpose of this research is to develop a simple but robust approach to 
modeling plate interactions using a contact problem in 2-D elastic finite element 
method. An individual plate is considered as a continuum, whereas an aggregat 
of plates is treated as a discontinuum such that plate boundaries are 
represented as contact surfaces. The behavior of plates at the contacts is 

defined hv the rmilnmh-X'avier failure criterion, and three types of contacts are 

considered: sliding contact(slip). tension release(separation). and sticking 

contact( single-node-continuum). The first and second modes correspond to 
transcurrent and divergent plate boundaries respectively. Convergent motion 
at plate boundaries is achieved by the double-node differential displacement 
technique with a recursive solution. Internal deformation caused by mid- plate 
processes, and body forces caused by elevation changes, are also implemented. 
A set of tvpical examples are discussed to validate this approach. To demon- 
strate its application, preliminary modeling of the Eastern Mediterranean, w ere 
tectonic deformation is produced by the interactions of the Eurasian, Arabian 
and African pLates, is carried out. An attempt is made to constrain the isp ace 
mpnt and stress fields bv geological observations and earthquake stress fields, 
” S pecUvely wl Kt the deformation pattern in the Eastern Mediterranean 
is substantially controlled by differential motion between the African and Ara 
bian plates and gravitational forces. Westerly escape of the Anatolian block 
accelerates markedly when buoyancy forces caused by elevation changes 
taken into account. 

Neotectonic Investigations (Appendices 2 and 3) 

Analysis of historical and instrumental earthquakes as well as recent field 
investigations of seismic surface fault offsets conducted by our group are usea 
to understand better the seismotectomcs of the North Anatolian fault, recent 
rates of plate deformation, the relationship between seismic slip and fault 
geometryf and the potential for future gap-filling earthquakes. These analyses 
suggest that recurrence intervals for large earthquakes is controlled b> e 
geometry and length of individual segments. Furthermore, the . east f e ™ p ^° f 
the westward/escaping Anatolian block appears to be divided in o g 

shaped blocks each of which moves independently. Recurrence intervals from 
historical earthquakes and geological data indicate a slip rate of 0.8-1 cm/yr or 
the eastern segment of the North Anatolian fault, suggesting abou - m sip 
deficit on the 1794 rupture segment. Thus, this segment is identified as a poten- 
tial seismic gap. 


GPS Measurements in Eastern Turkey (Appendix 4) 

During the Summer and Fail of 1989, MIT In cooperation with the Turkish 
Union of Geodesy and Geophysics (TUJJB), IFAG, and Durham University under- 
took a GPS field campaign aimed at establishing first epoch relative point posi- 
tloning in the region of continental collision in Eastern Turkey. As p&rt 
campaign, SLR sites were observed with GPS, footprints were established around 
each of the Turkish SLR stations, 16 regionally distributed sites were observed in 
Eastern Turkey, and a dense network was monitored in the Aegean Trough 
region of Western Turkey. Ten of the stations observed in 1989 were reobserva- 
tions of sites established and observed by MIT in 1988. Our group, in cooperation 
with our collaborators, plans a reoccupation in Western Turkey in 1990 and in 
Eastern Turkey in 1991. 


APPENDIX 1 


Finite Element Modeling of Plate Motions and its Application to the 

Eastern Mediterranean 


ORAL, M.B., Earth Resources Laboratory, Department of Earth, Atmospheric and Planetary 
Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA. 


Abstract 

Plate deformation patterns present a complex picture, and axe primarily affected by the inter- 
actions of plates at their boundaries and by their internal processes. The purpose of this study 
is to present a simple but robust approach to modeling plate interactions using a contact prob- 
lem in 2-D elastic finite element method. An individual plate is considered as a continuum, 
whereas an aggregate of plates is treated as a discontinuum such that plate boundaries are 
represented as contact surfaces. The behavior of the plates at the contacts is defined by the 
Coulomb-Navier failure criterion, and three types of contacts are considered: the sLiding con- 
tact(slip), the tension release(separation), and the sticking contact(single-node-concinuum). 
The first and second modes correspond to transcurrent and divergent plate boundaries, respec- 
tively. The convergent motion at the boundaries is achieved by the double-node differential 
displacement technique with a recursive solution. Internal deformation caused by mid -plate 
processes, and body forces caused by elevation changes, are also implemented. A set of typical 
examples are discussed to validate this approach. To demonstrate its application, a prelimi- 
nary modeling of the Eastern Mediterranean, where tectonic deformation is produced by the 
interactions of the Eurasian, Arabian and African plates, is carried out. An attempt is made 
to constrain the displacement and stress fields by the geological observations and earthquake 
focal mechanism solutions, respectively. The deformation pattern in the Eastern Mediter- 
ranean is substantially controlled by the differential motion between the African and Arabian 
plates and the gravitational forces. However, the westerly escape of the Turkish block accel- 
erates markedly when the buoyancy forces caused by elevation changes are taken into account. 
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Introduction 


Boundary forces/displacements caused by the interaction of plates at their boundaries, as well 
as internal deformation at mid-plates, gravitational forces and basal tractions, create extreme 
deformation. The finite element methods serve as a very applicable tool in understanding de- 
formational behavior of interacting plates. In analysis of plate deformation patterns, the 
geometry of the plate boundaries requires a certain attention: the creation and destruction of 
material, as well as relative plate motions are observed at the boundaries. Hence, this study is 
primarily engaged in for searching and investigating finite element models, and implementing 
the 2-D contact problem for moving plates, that incorporates pre-existing plate boundaries 
into solutions, in order to analyze deformational patterns caused by the interactions of the 
Eurasian, African and Arabian plates, and Turkish Block, where boundary displacements and 
gravitational forces are very significant. 


Since the classic paper by Turner, Clough, Martin and Topp which appeared in 1956, the 
deformation pattern of moving plates(bodies, mechanical parts, etc.) in contact, using finite 
element technique, has been studied by several authors in a variety of disciplines, e.g., earth 
sciences, civil and mechanical engineering, aerodynamics and geomechanics. Efforts in model- 
ing are also diversified in their approaches to the problem as well as their solution techniques: 
strong emphasis on the rheological behavior of the plates, concern for the behavior and for- 
mulation of the contact surfaces of the plates, and derivation of the governing mathematical 
relationships and their solution techniques. 

England and McKenzie(1982, 1983) suggested a thin viscous sheet model for continental de- 
formation which led to a number of studies (England, Houseman and Sonder, 1985; England 
and Houseman, 1986; Houseman and England, 1986; Sonder, England and Houseman, 1986) 
examining the deformational pattern of collision. They regarded the aggregate of plates as 
a continuum, which obeys a Newtonian or a power-law rheology. This approach ignored 
faulting/failure between the plates. Papers following this approach continued to define the 
medium as a continuum, and motions along the boundaries(e.g. faults), alternatively at the 
contacts of plates, are regarded as accommodating the strain rate field. 
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The crustal deformation in southern California, modeled as a creeping flow in a non-linear 
continuum(Bird and Piper, 1980) formed the basis of a discontinuity approach(Bird and 
Baumgardner, 1984). Both studies utilized the finite element method for solving their gov- 
erning equations. The former study assumed that the aggregate of plates is in a state of 
membrane stress subjected to plate- tec tonics boundary conditions. The flow law of this 
membrane contained a rigid-plastic term to represent frictional faulting in the upper crust 
and a power-law term to represent the dislocation creep in the lower crust. The continuum 
approximation to the region in which the strain-rate is fixed precludes the prediction of slip- 
rates on faults(contact surfaces). In the latter paper, the zones of contact(faults) are modeled 
by “special” elements to account for slip rates. 

In the engineering disciplines, emphasis on the behavior at contact first occurred in devel- 
opment of the joint element( 1968-Element) of Goodman et al.(1968). For the analysis of 
foundations and joint systems, Wilson(1976) employed an interface element. The analysis of 
jointed rocks(Goodman, 1976) was further refined to account for several modes of behavior 
at the con tacts( Goodman, 1975; Goodman and St. John, 1976). Wang and Voight(1969) 
alternatively utilized a contact algorithm to account for the behavior of the moving plates 
in contact. Goodman’s joint element method, using parabolic failure, defines the following 
modes at a contact: closing, opening, rotation, sliding with/without dilatancy/contractancy. 

In the finite element contact algorithm, used by Wang and Voight(1969) obeying the Coulomb- 
Navier failure criterion, which accounts for the behavior of the moving plates, three modes 
of behavior at a contact are defined: slip, separation and single-node-continuum. The state 
of stress and stability in underground openings is investigated using this contact algorithm. 
This approach is then applied to the progressive failure of rocks subjected to shear deforma- 
tion(Kasapoglu, 1973), to the analysis of contact surface traction due to glacial arcuate abra- 
sion cracks(Johnson, 1975), and to the collision of the Arabian and Eurasian plates(Kasapoglu 
and Toksoz, 1983). Further developments of contact algorithms are applications to 3-D static 
and dynamic analyses (Bathe and Chaudhry, 1985; Chaudhry and Bathe, 1986; Bathe and 
Mijailovich, 1987). They define the following modes of behavior at a contact: sticking and 
sliding contacts, and tension release. In a broad sense, the results of joint element and contact 
algorithm are equivalent (Bathe, 1989; pers. comm.). A hybrid finite element formulation 
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of the contact problem is achieved by Kumobora(1979) to investigate microscopic, moderate 
and extensive sliding. 

Following Wang and Voigt(1969), the contact problem is developed to explain all types of 
plate motions and to include gravitational forces and internal deformation. Also, a potential 
energy derivation of the governing finite element equations for this problem is introduced(see 
Appendix: Finite Element Formulation). Major emphasis is given to verifying the localiza- 
tion of failure zones, viz., slip-rates on the transcurrent plate boundaries(contacts). To avoid 
rheological complexities and disputes, and eliminate complicated mathematical treatise, a 
linear stress-strain relation for elastic media is assumed throughout the study. The following 
types of plate boundaries and behavioral modes on the contact surfaces are considered: 

1. For continent-continent collision, the contact surface is in sticking contact mode(single- 
node-continuum); 

2. For divergent boundaries, the contact surface is in tension release mode(separation); 

3. For transcurrent boundaries, the contact surface is in sliding contact mode; 

4. For subduction, the contact surface is in sliding contact(slip). 


Ideally, modeling of plate motions requires a 3-D analysis which is very voluminous and time 
consuming. However, introducing some a priori assumptions about the dip angle of the con- 
tact surfaces may pave the way for 2-D analyses. Assume that the boundary has a 90° dip 
for the first three cases and a 0° dip for the last one. Also, assume that material properties 
throughout the elastic plate are constant. Thus, dependency of the geometry in z-axis is re- 
movable. This consequently leads to 2-D analysis (Fig.l). To account for the sliding contact 
behavior along the subduction surface, a fourth mode( overlap mode) is defined. The calcula- 
tions are carried out for plane-strain case, using 4-node quadrilateral serendipidity elements 
with bilinear shape functions(see Appendix). Throughout the study, the term “displacement” 
is equivalent to instantaneous plate velocity. 

The input data, which consists of nodal point coordinates and element connectivity, is gen- 
erated by a pre-processor developed for this study. The results are displayed using a post- 
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processor(see Appendix). In the solutions of finite element equations, a versatile program, 
which includes the contact algorithm, is written. 

The next section gives a finite element treatise on continuum vs discontinuum approach. 
(Extensive treatise on the derivation of governing finite element equations are presented in 
Appendix.) After introducing contact algorithm and applying to a set of examples, the finite 
element models of present-day tectonics of the Eastern Mediterranean are discussed. 

Finite Element Models: Continuum and Discontinuum 

The continuum is defined as follows : 

Consider an aggregate of plates with various conditions imposed interacting at their bound- 
aries . If the pre-existing boundaries are not included in computations (no a priori boundaries), 
but inferred from the displacement /strain/stress pattern changes , which are spread over a re- 
gion rather than being localized , then this aggregate is considered as an individual plate which 
deforms under the combined effects of various conditions posed on each member. 

To demonstrate this, interactions of two plates are considered. Their total area is 700 km x 
500 km. Their boundary which makes an azimuth with the north. In Fig.2, the deformation 
pattern for the simulation of a transcurrent boundary is given. Only northward edge dis- 
placements (instantaneous plate velocities), 1 cm and 3 cm relative to a reference plate, are 
applied to the lower boundary of plates A and B, respectively (Fig.2a). To avoid any struc- 
tural instability, rollers are applied to prohibit any east- west motion. This also accounts for 
the sideways continuity of the plates. The expected transcurrent plate boundary is denoted 
by dashed lines and meshed with smaller elements. Note that a definition of two distinct 
plates is made by a sudden change in boundary conditions. The occurrence of slip along the 
boundary so that one plate slips past the other is anticipated. The deformed structure plot 
(Fig.2b) shows the motion of the plate(s). However, a close look at the displacement vector 
configuration(Fig.2c), shows that the anticipated failure is distributed over a wide range in 
the east-west direction rather than being localized and that there is no localization. It is 
obvious that the differential displacement(slip) is linearly distributed between the terminal 
ends of the plates. Though the structure is bounded on its lateral sides and compressed on 
the lower side, tensional features develop. There is no observable stress difference ( drop) 
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in the vicinity of the boundary. The maximum shear stress fade out is almost identical on 
both sides. NE-SW expansion and NW-SE contraction is observed. The strain around the 
boundary is markedly large and it decays with distance. 

Fig.3 simulates of a divergent plate boundary. Using the mesh for the previous model, pulls 
of -1 cm and 1 cm are applied to east and west boundaries of the plates A and B, respectively 
The deformed structure plot and displacement pattern, shows no opening, but stretching in 
the E-W direction and a small amount of shortening in N-S direction. Stress and strain pat- 
terns give no clue for the possible location of failure. 


Next, these two cases of deformation are treated as a contact problem and solved in discon- 
tinuum with the contact algorithm (Figs.3 and 4). Details of the contact problem will be 
discussed in the following section. The discontinuum is defined as follows : 


When an aggregate of plates with various conditions imposed on them interact at their bound- 
aries such that the pre-existing boundaries are included in computations using dual nodes 
(not split nodes) , and their behavior at the contacts is determined by a frictional law , this ag- 
gregate establishes the discontinuum , and the individual plates are considered as a continuum 
deforming under the conditions prescribed . 

The boundary conditions for the example in Fig. 2 are used for a similar mesh in which the 
boundary between two plates is shown by a solid line and the deformation field is solved as 
a contact problem. In Fig.6c (cf. Fig.2) there is a marked slip (~ 2cm), which is the differ- 
ence between the motions of the two plates along the pre-defined plate boundary (Fig.6a). 
This illustrates sliding-contact behavior. There is a remarkable shear stress drop (Fig.6e) 
across the fault. Because no motion is allowed in the E-W direction, compression develops 
perpendicular to the fault. The NW-SE compression implies that the differential motion is 
translated across the boundary because the slower plate has nowhere to go but north. As a 
result, principle major and minor stresses in plate A, in comparison to plate B increase and 
decrease, respectively. The slower plate (A) stretches almost sub-parallel to the boundary 
while it shortens perpendicular to it. Note that the stress and strain patterns change with 
the length scale of the plates in the E-W direction which are controlled by the strike of the 
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boundary. 


For a comparison, the divergent plate boundary case is treated in discontinuum as a contact 
problem. In Figs. 5b and c (cf. Fig.3) the opening between the plates is 2cm which is equal 
to the differential motion between the plates. This demonstrates the tension-release-contact 
behavior. Strain for both of the plates is zero as is the stress field, since the displacement 
field is constant. 


The magnitudes of strain and stress will be discussed before other finite element models are 
analyzed. 

Contact Problem 

The problem of pressure distribution between two bodies was first solved by Hertz(1881). In 
the Hertz elastic theory of contact (Timoshenko, 1934), the contact surfaces are frictionless 
and do not transmit tangential surface tractions across their boundary. However, Cater(1926) 
showed the inevitable occurrence of sliding within the area of contact. This study uses fric- 
tional law defined by the Coulomb- Navier criterion. ( Its incorporation into contact problem 
will be defined later.) A contact happens when two or more bodies meet each other at a con- 
tact surface (Fig. 6) is created. Impact is generated when a dynamic contact occurs (Johnson, 
1976). 

The forces which develop at the contact surface determine the behavioral mode. Unless there 
is a tensional release, the forces acting at the contact must be equal and opposite in sign, 
otherwise they are zero. Moreover, the relation between tangential and normal forces acting 
on the contact surface must satisfy certain conditions determined by the Coulomb- Navier 
criterion (Jaeger and Cook, 1979). Accordingly, failure occurs when the shear stress exceeds 
a threshold defined by normal stress scaled by a, factor and cohesion of the material : 

II T ||= S 0 (1) 

where r, <7, S 0 and p are the shear and normal stresses, shear strength(cohesion) and friction 
coefficient (Fig. 7). Therefore, slip will occur when 
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and tensional release will take place when 


(7 > T 0 , (3) 

where T 0 is tensile strength and as defined positive. 

The contact forces developed on the contact surface can be related to cr and r by averaging 
them over the distance C half-way in either direction from the contacting node which lies on 
the contacting segments(Fig.8): 


if” (4) 

and 

§=IM|, (5) 

where P n and P t are normal and tangential contact forces and their relation to P x and P y is: 

P n = P x cosa + P y sina 

( 6 ) 

Pt = — P x sina + P y cosa, 

where a is the angle between the normal to the contact surface and the x-axis (Fig.8). In 
finite element mesh and computations, the boundary(contact surface) of contacting plates is 
represented by dual nodes, and the forces which develop on this surface are taken into con 
sideration when judging the behavioral mode at the contact. If the contact forces are totally 
ignored, dual nodes become split nodes. This is a violation of compatibility in finite element. 
Four modes are defined: sticking contact, sliding contact, tension release and overlap. Mode 
IV is introduced to account for subduction and will be discussed later. Analyses of the first 
three modes are given below. 


Mode I : Sticking Contact(single-node-continuum) 


When contact stresses are insufficient to activate any motion along the contact surface such 
that 


Pn <T 0 C, 

Pt < So C - fl Pn , 


then the components of contact forces satisfy 


(7) 


P? + P* = 0 , 
Pf + P* = 0 , 


( 8 ) 


where superscripts denote the plate at the contact surface. Hence, the displacements are 




= 


= ?? 


(9) 


so that the contacting plates in discontinuum behave as if they are in continuum. 


Mode II : Sliding Contact(slip) 


The plates in contact exhibit a sliding behavior when tangential stresses exceed the linearly 
defined Mohr envelope, provided that the normal contact forces are less then the tensile 
strength: 


Pn <T 0 £ , 

P t > So C - \i P n . 


( 10 ) 


The amount of slip is determined by the frictional coefficient, and shear and tensile strength. 
From the continuity of the stresses across the element, the components of the contact forces 



with additional constraints by the frictional law: 


II P A II -So C + n P A = 0 , 

II A S II -S 0 C-nPB= 0. 

This states that the displacements perpendicular to the contact surface are equal: 


( 12 ) 


implying that 


9* + Qx = tana [?y + 1 ' 


(13) 


<7x = . 

9y =«?’ 


for a = 0° , 

for a = 90° 


(14) 


Mode III : Tension Release(separation) 


The plates diverge when the normal stress exceeds the tensile strength: 


Pn <T 0 C , 

in which the components of the contact forces are 


(15) 


Pf =0, 

Py A =0, 

Pf = o , 
Pf = 0 - 


(16) 


Note that (a) the displacements(at the nodes and contacts-dual nodes) are unknown, (b) the 
contact forces determine the behavioral mode at the contact, (c) the contact forces also are 
unknown. To settle the dust, recall Eqn.A14 and rewrite as 


Kq — Q = —P , 


(17) 
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where P contains P x and P y unprescribed contact nodes. The contact forces P can be calcu- 
lated and used to determine the behavioral mode, provided that the nodal displacements are 
known. This requires another iterative scheme in which an interrogation procedure (Wang 
and Voight, 1969) for behavioral modes is constructed. It is clear that P = 0 when no contact 
is defined and the displacements can be calculated from Eqn.A42 in conjunction with Eqn.A41. 


Calculation of Displacements at Contact Surfaces: A Contact Algorithm 


To achieve a good convergence, the initial displacement vectors have to be close to those of 
the final solution. A first approximation can be drawn from a sticking contact solution, where 
discontinuum becomes asymptotic to the continuum. For the first part of iterative solution 
of the displacements at dual nodes, the following equations prevail(see Appendix): 


q\ - q[ 1 - u {<p\ + + ^i+l.i+l) , 

<?«+! =?'• 

Displacements at the other nodes are obtained from Eqn.A42. 


(18) 


In the rest of the iterative solution, the contact algorithm is applied. To decide on the 
behavioral mode, the displacements at one of the dual nodes have to be calculated from the 
relation: 


Ki,i + Ki + 2,i+2 Kt,t+l + K i+2, {+3 

9» 


-(v»! + A+i) 

A'i+1,,- + Ki+Z'i+2 A'i+1,,+1 + +3,j +3 

. 9<+l . 


. 1 + ^'+ 3 ) . 


where q, and g,+i are x- and y-parallel displacements at dual node k on plate A. In / - f/i 
iteration, the contact forces equal the residual defined by Eqn.A41: 


P l n = <p l cosa -(- (A sina , 

Z I ■ /, ( 20 ) 

Pi = -<fix stna + <f y cosa , 

For the contact algorithm, the angle between the contact and its normal a, the distance that 
the nodal contact forces uniformly distributed C , and the trigonometric relationships are: 
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% 


c = vT ~k + **+ IF + (y* - yk+TF , 


fana = 
sina = 
cosa = 


Vk-Vk+i 

-Xk+Xk+l 

-Vk+yn+i 


y 


y 


y 


( 21 ) 


where x and y are the coordinates of the k — th dual node. Dual nodes, A and B may or may 
not share the same location. For the purposes of this study, they have identical coordinates. 


For mode I, no further calculations are necessary, and the displacements on either side of the 
contact are equal: 


</i+2 = ?! 

9i+ 3 = 9»+li 

where g,+ 2 and <7,4.3 are x- and y-parallel displacements at dual node k on plate B. 


( 22 ) 


Mode II calculations require rotations of diagonal elements of the stiffness matrix in the 
direction normal and parallel to the contact surface, to account for the tangential and normal 
stresses(see Eqn.l2b). Define 


xa — — sina ± (—ll) cosa , _ 

v ; (23) 

ya = cosa ± (— fi) sina , 

where xa and ya are rotation coefficients. If P t > 0 ” is used. To solve for q { , qi+u <7,4-2 

and qi+ 3, first set them to zero( which removes the force contribution at these nodes, recall 
Eqn.A41) and then solve : 


K{,ixa + Ki+\',ya 

Ki+\'i+iya + Ki'i+\xa 

0 

0 


?' 1 

K, ti 

Ki , i+ 1 

Ki+l'i+2 

^*+2,1+3 


9 <+ 1 

K i+U 


K i+3,i+2 

^t+3,1+3 


?!+ 2 

1 

tana 

-1 

tana 


. 9 , >3 . 


-(fact + <p\ +1 ya - S 0 C) 
~(fi + ^'+2) 

-(vi+i + ^+3) 

0 


(24) 
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For mode III, there is no bonding between the dual nodes lying on the contact surface, and 
the displacements for each plate follows a different linear system solution. To solve for g, and 
<?,■+! ; first set them to zero and solve: 


Ki,i 

*.,+1 




-v{ 

. Ki+U 

-^+i,i+i 


. ^'+1 . 


. “Vi+1 . 


And to solve for qi + 2 and <7,4.3; first set them to zero and solve: 


(25) 


A\+ 2,t+2 #»+2,i+3 

A\+3,,'+2 ^i+3,1+3 

The calculation of displacements for the rest of the nodes is carried out by Eqns.A42 and 
A41. Convergence in the first stage of iterative solution is monitored by Eqn.A41. In the 
next stage, 


tfi+2 


*»+3 


“Vt+2 

M+3 


(26) 


C = £ llvJ-ll. (27) 

» 

where M is the set of degrees of freedom involved contact. C is the measure of convergence 
and must decrease as more iterations are performed. 

On the Magnitudes of Strain and Stress, and the Strength and Coefficient of Friction 

The average strain per year is at the order of. 10 -7 (Turcotte and Schubert, 1982). Lithostatic 
stress for a 35 km thick plate is 1000 Mpa (Turcotte and Schubert, 1982). The compressive 
stresses caused by elevation changes(Molnar and Lyon-Caen, 1988), and the stresses induced 
by ridge-push and slab-pull (Richardson, 1972), range from 60 to 33-200 MPa, while thermal 
and membrane stresses are 400-600 Mpa. Bott(1982) reports an average of 100 Mpa. 

Using the instantaneous plate velocities, the magnitudes of the strain(rate), found in this 
study, are at the order of 3 x 10 -7 coheres with the reported values. Stress magnitudes are, 
however, at the order of 3000 Pa(= 12 GPa x 3il0 -7 ). Compared to the magnitudes given 
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above, there is a scaling factor of 10 3 - 10 4 which arises because (a) instantaneous plate ve- 
locities axe used as boundary displacements, (b) strains calculated from these displacements 
are strains per year(strain rate), and (c) stresses calculated from these strain rates inherit 
this time scaling and are stresses per year. 

This “time-factor” x also scales the contact forces which means that the magnitudes of shear 
and tensile strengths are also scaled corresponding ~ 1000 — 10000 years(episode cycle). The 
calculations using zero strength do not significantly differ from those using this value, when 
the forces are greater than the strength of the plate. Higher strength values cause slowing in 
plate motions. Unless otherwise stated, all models have zero strength. For the sliding mode 
this implies sliding without cohesion. 

The friction coefficient is another factor that reduces the motions along plate boundaries. 
Bird and Baumgardner (1984) point out the possibilities of a low friction coefficient (~ 0.3) 
for active faults. It is found that fi < 0.5 produce almost the same deformation pattern. 
Throughout the study, a “low” friction coefficient is used. 

No strength and a “low” friction coefficient are appreciable (a) when the contact surfaces 
have already formed and are weak, viz., have an extreme tendency to move, (b) when the 
instantaneous motions are modeled, assuming a certain amount motion along the boundaries 
per year, viz M “things change”. 

Examples for Plate Motions 

Contact solutions of interacting two plates are demonstrated (Figs. 4 and 5) and compared 
to those of continuum solutions. This section depicts some aspects of contact solutions for 
three-plate cases. At the contact of three plates, triple nodes are required to simulate the 
discontinuum. However, this may be avoided when the third plate boundary is translated 
one element away from this junction. The following models exploit this numerical treatment. 
Note that the goal of the study is to demonstrate the advantages of the contact problem, 
so a more general algorithm is not incorporated. When n-plates are in contact, instead of 
using contactor nodal point forces (this study; Wang and Voight, 1969), contactor segment 
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tractions (Bathe and Chaudhry, 1985) calculations become more attractive. 

In the three-plate-contact models, the following considerations are in effect: (a) a roller sub- 
parallel to N-S or E-W implies that there is no motion allowed in their normal directions, 
and this boundary corresponds to transforms; (b) a hinge means that there is no motion in 
either direction, which in turn serves as a reference frame; (c) there is a hypothetical ref- 
erence frame for the instantaneous plate velocities (displacements) located somewhere away 
from the plates unless otherwise located; (d) the contact algorithm previously discussed is 
exploited; (e) plates are named as lower (right) plate A , lower (left) plate B> and upper plate 
C; (f) at the end points of contacts, there is incompatibility in the deformed structure plots 
when either one of the components of dual nodes is constrained, or the requirement for triple 
node is necessary; (g) unless otherwise stated, there is no overlap, this being controlled by 
additional constraints(double-node differential displacements). This kind of overlap occurs 
due to the effects mentioned in (f) and since linear interpolation in joining two points is used, 
these effects are enhanced. For those nodal points, refer to displacement vector plots. 

In Fig.9, the interactions of two plates with a fixed reference plate is considered. On the 
boundaries of right and left plate, 3 cm and 1.4 cm displacements in the N and NW direc- 
tions are applied, respectively. Most of the differential motion is taken between the right and 
left, and between right and upper plates, as sliding contact. This type of model could be 
considered a Fault-Fault-Fault type triple junction. The decrease in slip along the right plate 
contact is mainly compensated by internal deformation as shown by the strain field (Fig.9f). 
The motion of the faster plate decreases toward the north for (a) there is no sideways motion 
allowed and (b) it meets another stable plate. This stress pattern demonstrates the sub- 
groupings of principal directions. Compression is dominant in the upper plate as well as in 
the right plate. The upper plate also has a tendency to extend sub-parallel to the contact 
surface with the lower right plate, due to the loading of the right plate. A striking change in 
maximum shear stress magnitudes marks the contacts. 

For the model shown in Fig. 10, the boundary conditions are identical to the previous one 
except that the upper plate is no more fixed, but is displaced towards N by 1 cm. These 
boundaries could be considered as Fault-Fault-Ridge triple junction(Fig. 10b). The differen- 
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tial motion between left and upper plates is not enough to create the tension release mode. 
However, their interaction, especially the motion of upper plates, triggers divergent plate mo- 
tion. Compared to Fig. 9, the slip along the contact between upper and right plates decreases, 
whereas it remains the same along the contact between right and left plates. Due to ten- 
sion release, the northerly oriented stresses in the left plate vanish while the E-W extension 
increases. Evidently, the change in maximum shear stress patterns occurs at plate boundaries. 

The deformation pattern in Fig. 11 is an example of modes I and II, sticking and sliding 
contacts, respectively. Boundary conditions are the same as in the previous model ex- 
cept that the left plate is displaced northernly by 1 cm and set free on its west end, and 
that the sense of upper plate motion is reversed. It may be considered as a Fault-Fault- 
Trench(collision) triple junction. The contact between upper and left plates is the first type 
of plate boundaries(continent-continent). The displacement field seems to be zero at their 
contact because of the opposite polarity of the plate motion. One of the dual nodes at the 
east end is not constrained and causes an artifact motion, and has to be ignored. The entire 
region is in compression. When stress and especially strain fields are considered, the left and 
the western of the upper plate behaves as a single plate in continuum (sticking contact). As 
a corollary to this observation, the maximum shear stress pattern gives no information about 
the existence of a contact. 

The boundary conditions for the model in Fig.12 are the same as in Fig.10, except that the 
left plate is displaced northwesterly by 1.4 cm. To simulate Fault-Fault-Trench(subduction) 
triple junction, mode IV is introduced: 

Overlap mode 


Consider nodes q* and qf and let 6 be the convergence amount between these plates. 
6 is the measure of material subducted and could be correlated to dissipation of mate- 
rial(displacement) over a region such that: 


Qs 


= < 7 ^ + ? 


(28) 
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where qf and qf are displacements on plates A and B at nodes s and r, respectively and 
related to each other by a constant of dissipation(rate of subduction) 6 . A first solution may 
be obtained by ignoring convergence, and then incorporating the above relation, an additional 
set of prescribed degrees of freedom may be obtained. To simulate this type of plate boundary, 
a recursive solution could be using the additional constraints and former boundary conditions. 
This solution is called double-node-differential displacements with recursive solution(Overlap, 
Mode IVa). 

A tentative solution is given in Fig. 12 with boundary conditions the same as Fig. 10 except 
that the sense of motion in the upper plate is reversed. 3 mm convergence is included across 
the contact of the upper and left plates. The displacement field demonstrates the success of 
the approach. Because of this motion, the principal stresses normal to the overlapping contact 
drastically decrease. The change in maximum shear stress pattern corresponds to the overlap. 

An alternative solution can be obtained by considering the auxiliary contact, perpendicular 
to the primary contact, and then applying double-node differential displacements( Overlap, 
Mode IVb). This removes locking, and perfectly transmits the motions from one plate to the 
other. To avoid the very complicated mesh required for an auxiliary contact, the solutions for 
convergent motions, at this stage , are carried out by double-node differential displacements 
using primary contact. 

Internal deformation caused by folding, kinking, etc., is modeled and shown in Fig. 13, where 
the boundary conditions are the same as the previous example and exclude convergence. The 
elements with crosses are assigned negative initial strain corresponding to (1cm) shortening. 
When this region becomes softer, it takes up some of the differential motion between the 
plates. The displacement pattern shows that the slip along the boundary between the left 
and right plates decreases. This contribution boosts the stress and strain field over this re- 
gion. The same result is obtainable by decreasing the Young modulus for these elements. 

McKenzie(1978) discusses three possible driving forces: boundary forces which are widely 
used until now in the above examples; gravitational ( buoyancy )forces, caused by elevation 
changes; and forces on the base of the lithosphere. When the area of the plates is rela- 
tively small compared to boundary and gravitational forces, basal tractions may be ignored. 
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Therefore, they are excluded in finite element computations. However, body forces caused by 
elevation changes are very significant, and included. Molnar and Lyon-Caen(1988) calculated 
the upper bounds these buoyancy forces. 

Fig. 14 shows the effect of body forces on the deformation pattern. The boundary conditions 
are identical to those in Fig. 11 except that the east end of upper plate is free to move in the 
E-W direction. The left and upper plates are combined and allowed to behave as a single 
plate. Although this contact is locked, the upper plate escapes westwardly through the free 
end due to the applied body force. The exact amount of body forces is inconclusive since no 
lower bound is defined. In this study the first estimates are obtained from the Molnar and 
Lyon-Caen’s calculations for horizontal driving forces. After normalizing per unit volume and 
per unit year, a scaling factor of 10-100 is observed. 

Application to the Eastern Mediterranean 

The region of interest includes the Eurasian, Arabian and African Plates, as well as the 
Turkish and East Anatolian Blocks. The neotectonics of the region is shaped by the second 
opening episode of the Red Sea during the early Pliocene (4.5 Ma). The Arabian Plate, for- 
merly having the same velocity as the African Plate, gained acceleration and, following the 
closure of the Thetys Ocean, collided with the Eurasian Plate. The African Plate, on the 
other hand, continues its subduction under the Hellenic and Cyprean Arcs. This continuing 
continent-continent collision between Eurasia and Arabia creates extreme deformation in the 
region. Shortening in eastern Turkey was first accommodated by crustal thickening. Later, 
instead of excessive crustal thickening, the Turkish Block which is bounded by the North and 
East Anatolian Faults, wedged out towards the west under the compressive regime. Presently, 
the western tip of this block extends in a N-S direction, by accommodating the area of the 
African Plate lost by subduction(Fig.lS) as westerly motion is inhibited by the Grecian Shear 
Zone. The major structures in this region are the Dead Sea, the East, North and Northeast 
Anatolian Faults, the Cyprean Arc and the Bitlis Suture. At the Maras and Karliova triple 
junctions the Arabian, African and Turkish, and the Arabian, Turkish and Eurasian plates 
meet. Fig. 16 summarizes the sense of motion along the faults. Accordingly, the sense of 
displacement on the Dead Sea, and the East and Northeast Anatolian Faults is sinistra!, 
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whereas on the North Anatolian Fault it is dextral. Several authors calculated the rota- 
tion poles of these plates (Chase, 1978; McKenzie, 1972; Minster et al., 1978; Gordon and 
Jurdy, 1986). From these data, local plate velocities(Fig.l7) are computed for the African 
and Arabian plates with respect to the Eurasian plate(Cox and Hart, 1986). Local plate 
velocities obtained from McKenzie(1972) are higher than the others. Mean boundary dis- 
placements(instantaneous plate velocities) for the African and Arabian Plates are northerly 
5.2 mm and 22 mm, respectively. Fig. 18 shows the area where finite element calculations are 
applied to the Eastern Mediterranean. In the west, it is bounded by the Pliny-Strabo Trench 
system, and excludes western Turkey. The Dead Sea, East, North, and Northeast Anatolian 
Faults, the Cyprean Arc and the Bitlis Suture are considered and modeled as 2-D contact 
surfaces. 

Three finite element models of the Eastern Mediterranean are investigated. They all consider: 

(a) strike-slip faulting[model l](Fig.l9), 

(b) strike-slip faulting and convergence at the Cyprean Arc and Bitlis Suturefmodel 2](Fig.20), 
and 

(c) strike-slip faulting and convergence at the Cyprean Arc and Bitlis Suture, internal defor- 
mation at the Palmyra Kink and gravitational forces that wedge out the Turkish and East 
Anatolian Blocks[model 3](Fig.21), 

to be major tectonic elements in creating the deformational pattern. The first model has only 
historical importance and proves the need for the latter two models. At the beginning, the 
finite element models of the Eastern Mediterranean only accounted for strike-slip faulting and 
ignored convergence, especially the subdue tion of the African plate at the Hellenic Trench 
and the Cyprean Arc(Kasapoglu and Toksoz, 1983; 1988). These models, unfortunately, cre- 
ated skepticism about 2-D finite element calculations in a region where some plate motions, 
generated by the differential motion between the African and Arabian Plates, is taken up 
by convergence(subduction/collision). Strong emphasis on the contribution of gravitational 
forces also invited speculations. Other than insisting suggestions to incorporate these forces 
in finite element calculations, no models have appeared yet. To settle the dust, convergence 
and body forces, as well as internal deformation, are included into 2-D finite element compu- 
tations. The contact problem previously discussed is utilized. A standard mesh is generated 
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using a pre-processor program developed for contact problems. The boundary conditions are 
identical in each model. Rollers at the sides of the plates illustrate the direction in which the 
relevant nodes may move. Rollers west of the African and east of the Arabian Plates reflect 
the uniform continuity of the plate motions beyond the modeled area. The Eurasian plate is 
held fixed. 5.2 mm and 22 mm northerly boundary displacements(instantaneous plate veloc- 
ities relative to the Eurasian Plate) are applied on the lower end of the African and Arabian 
Plates, respectively. The tension release mode is applied when one or both of the dual nodes 
at the terminal ends of plate boundaries are constrained causing incompatibility at these 
nodes. Therefore, the displacement/stress/strains at these nodes are ignored. In the second 
and third models, 3 mm convergence across the Cyprean Arc and 6 mm shortening at the 
Bitlis Suture is included. The third model includes 4 mm internal deformation taken up by 
Palmyra Kink(Barka et al., 1989) and gravitational forces. The gravitational force per unit 
volume is estimated from the upper bounds of the horizontal driving force(Molnar and Lyon- 
Caen, 1988). Assuming N-S direction components are balanced, westerly and easterly forces 
are applied on the elements that represent the Turkish and East Anatolian blocks, respec- 
tively(Fig.21a). A comparison of displacement fields of models 1 and 2 show that introduction 
of convergence at the Cyprean Arc and shortening at the Bitlis Suture, results in decrease 
of slip along the fault zones. Because these two features take up ~ 3 — 5mm, the motions 
along the North, Northeast and East Anatolian faults are halved. The westerly and easterly 
escapes almost vanish. The slip along the Dead Sea Fault decreases northerly in both models. 
These data imply that the differential motion between the African and Arabian Plates are not 
sufficient to explain the amount of slip observed along the Anatolian transforms. Barka and 
Gulen(1989) argue that the slip along the North Anatolian Fault is 1 cm/yr, if the age of the 
fault is early-middle Pliocene(3.5-4 Ma) and has a total offset of 35 km, which decreases west- 
erly to 5 mm/yr. The East Anatolian fault has a 5 mm/yr slip(Barka et al, 1989). Geological 
studies on the Dead Sea Fault show that the 10-15 mm/yr slip in the south(Gharb Segment) 
reduces to 5 mm in the north(Karasu segment). Given this information, finite element mod- 
eling becomes inconsistent with the observational data. “To do justice to the whole data 
set”(Sengor, written comm.), gravitational forces and internal deformation at Palmyra Kink, 
as well as convergence, are included into the computations( model 3). Fig.21 demonstrates 
how the Turkish and East Anatolian Blocks wedge out. The openings at the KarLiova and 
Maras triple junctions are reflected in displacement and stress patterns. On the average, 15 
mm and 5 mm slip along the Gharb and Karasu segments of the Dead Sea Fault are observed. 
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As some part of the differential motion is accommodated by the Palmyra Kink, build up of 
strain/stress occurs. The average slip rates on the North(13 mm in the east, 8 mm in the 
west) and Northeast Anatolian Faults(4 mm) are relatively large requiring a smaller body 
force magnitude. The primary conclusion is that the slip rates along the Anatolian trans- 
forms are controlled by the buoyancy forces caused by the differential stresses between higher 
and lower elevations. The gravitational force used for these calculations is lower than the 
upper bound for 2000 m elevation, but the largest for the admissable deformation patterns. 
When higher values are preferred, the tectonic picture changes: The transcurrent motion 
character along the East and North Anatolian Faults’ changes to divergent and convergent 
type boundaries, respectively. Thus, slip rates along the North Anatolian Fault is no greater 
than 1 cm/yr. Major stress concentrations take place in the Turkish Block, Palmyra Kink 
and north of the Bitlis Suture. Compressive stresses are characteristic to the Cyprean Arc 
and the Bitlis Suture. Because western Turkey and the Aegean Sea are not included in the 
models, transition to extensional regime is not observable. However, the westerly escape of 
the Turkish Block is reflected as a SW-NE compression. Strain patterns suggest that the 
strain accumulates north of convergent zones and is at the order of 1-5 x 10 -8 . 


Conclusions and Discussions 


The deformation pattern of interacting plates can be modeled using 2-D finite element method. 
The comparisons between continuum and discontinuum approaches demonstrate that the 
contact algorithm is a robust method for modeling plate motions. Derived potential energy 
expressions show that the contact problem is inherent in finite element equations and that 
implementing rheologies, other than linear elastic type, is straightforward. Application to 
the Eastern Mediterranean, where extreme deformation is created by the interactions of the 
Eurasian, African and Arabian Plates, shows that the regional tectonic picture cannot only 
be defined by boundary displacements(the ridge push force due to the opening of the Red 
Sea). Including gravitational forces is a must. These two forces control the deformational pro- 
cess in the region, and suggest that the differential motion between the Arabian and African 
Plates is responsible for the slip rate along the Dead Sea Fault, whereas the gravitational 
forces create the slips along the Anatolian Transforms. The models indicate that the slip 
rate must not exceed 1 cm/yr along the North Anatolian Fault. Accounting for the Palmyra 
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Kink demonstrated that mid-plate internal deformations are another source of strain accu- 
mulation. Introduction of the overlap mode shows that unless subduction and shortening are 
taken into account, understanding the deformational behavior of plates is incomplete, and the 
patterns obtained may lead to false conclusions. The overlap mode IVa illicits considerable 
success, however, better results can be obtained by using mode IVb. Further developments 
of the contact algorithm, viz., extending to the third dimension, triple-node formalism, us- 
ing contactor segment tractions in calculation of contact forces instead of nodal point forces, 
and employing the joint element approach; and the incorporation of other rheologies into 
the computations will improve modeling of the deformational behavior of (micro)plates with 
complicated boundaries. Instead of the assumed displacement finite element approach, it is 
suggested that assumed strain method, for regions where strain rates are significant and hy- 
brid finite element method for regions where prescription of both stresses and displacements 
are needed, be utilized. As long as the instantaneous plate velocities are used , which in- 
evitably introduces a time scaling, the low coefficient of friction and no-strength is justifiable. 
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APPENDIX : Finite Element Formulation 


Under boundary and nodal displacements, body forces, initial strain and stress and surface 
tractions, the deformation pattern of the plates is defined by displacement, strain and stress fields. 
Since instantaneous plate motions are considered in this study, they, hereinafter, will be referred 
to as displacements. A plane strain/stress media with linear stress-strain relation is assumed in 
the general derivation of governing finite element equations, including the contact problem. First, 
the finite element formalism using total potential energy is introduced. After finding the governing 
equations by variational calculus, discretization of the displacement field is discussed. Following 
stress calculations from strain field, numerical solutions techniques are given. 

Governing Finite Element Equations 

The linear relation between stress and strain for elastic materials in 2-D is given by: 

a = Ee - Ee 0 + a 0 , (1) 

with 


O — [ <T X (Jy 0 X y ] > € — [ € x € X y ] , 


( 2 ) 


where tr, e c„ and a 0 are stress, strain, initial strain and initial stress, respectively. For plane stress 
and strain, the material property matrices E are: 


E = 


Z 

1 - i/ 2 


1 v 0 
v 1 0 

0 0 ^ 


E = 


Z 

(1 - 2i/)(l + v) 


l-i/ 1 0 

0 1 - v 0 

0 0 


(3) 


respectively, where Z and v are Young’s modulus(~ lOOGPa) and Poisson’s ration (~ 0.25) strain. 
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Consider a conservative system on which external and internal forces act. The internal work 
done(the stored energy) U 0 of this system, per unit volume, is: 

Uo =f<Td(, 

= \i T Et - ( T E( + ( T <j 0 , 
and the total strain energy U is given by: 

U = J v U 0 dV + q T P, (5) 

Introducing the contact potential, the second term accounts for the contribution of the contact 
forces. The total external work done on this system is: 

W = J u T F dV + J u T <j> dS + q T R , (6) 

where P , u(u(x, y), v(x, j/)), F , <£, q and R are internal nodal forces(discrete)[contact forces], dis- 
placement (continuous), body forces, surface tractions, displacement(discretized) and concentrated 
loads(nodal point forces), respectively. Then, the total potential energy functional is 

n p = u-w. ( 7 ) 

The above equation indicates that the total potential energy is not simply a function of displace- 
ments and their derivatives, but also depends on their integrated effect. 


For a linear elastic medium, the strain can be related to displacement by a differential operator 
D , 


e = Du , (8) 

and the displacement field is discretized at the nodal points of an element by shape (interpolation) 
function matrix N, 

u = Nq , (9) 

and finally, an expression for strain is: 

e = Bq, ( 10 ) 
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where B is the strain-displacement matrix and has k~ 1 units when displacement q has k units. 
(The explicit expressions for D y N and B will be given later.) This procedure invokes assumed- 
displacement finite element method where the displacements at the nodal points represent the 
degrees of freedom (unknowns). In this method, the essential boundary conditions are prescription 
of displacements, and the natural boundary conditions are prescription of stress. This implies that 
prescribing displacements at the boundaries alone is sufficient to solve the equations. 

Substituting (4), (5) and (9) into (6), the toted potential energy becomes: 

n p = \q T (Z f v B T EBdV)q - «z r (£ f v B T Ee 0 dV ) 

-? T (E fv N T a 0 dV) - q T (Z f v N T F b dV ) 

~q T (Z Is N T <fidS) - q T R + q T P (11) 

= \q T (H k)q - q T {‘ 1 , 0 + 7<To + 76 + 7» + R - P) 

= \q T Kq - q T (Q - P) . 


Summation has to be done over the total number of elements indicating the requirement for the 
assembly of element stiffness matrices and nodal forces, k is the element stiffness matrix whereas 
K is the assembled globed stiffness matrix. Q represents the nodal forces caused by initial strain, 
initial stress, distributed forces and surface tractions. The stability of the equilibrium state is 
realized when the potential energy is minimum, i.e., when it is stationary with respect to u small 
variations” of displacement: 


lip = stationary . (12) 

This is the principle of minimum potential energy. According to the calculus of variations, the first 
variation of total potential energy with respect to q must be zero: 

£II p = 0, (13) 

such that: 

Kq = Q-P. (14) 
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This result is equivalent to that obtained by using the principle of virtual displacements(Zienkiewicz, 
1986) and is the governing equation for contact problems. The structure stiffness matrix K relates 
nodal point displacements to nodal point forces. P is zero for continuum where there is no contact 
defined; therefore, for continuum the governing finite element equation is: 

Kq = Q. (15) 

Thus, it is implied that Eqn.14 is the governing equation for discontinuum. Regardless of rhe- 
ology, the finite element equations would be reduced to one of these equations, depending on what 
kind of medium(continuous/discontinuous) is chosen. It is, therefore, inherent in the finite element 
approach to solve for the deformation pattern of discontinuous media. 

The j-th column of the stiffness matrix is the vector of nodal forces applied to maintain static 
equilibrium when j-th degree of freedom has unit displacement and others have zero displacement. 
The diagonal elements of the stiffness matrix are positive. A zero diagonal element would create 
zero reaction force Q, creating an unstable structure. It cannot be negative because this would lead 
to a physically unrealizable situation where displacement and force vectors lie in opposite direction. 
If only linear degrees of freedom are considered, the sum of each column is zero. When there is 
a linear relation between force and displacement, the stiffness matrix is symmetric. The stiffness 
matrix, however, is singular. Its rank is less than its dimensions by the number of rigid body 
motions. Eigenvalue analysis of the stiffness matrix determines the number of rigid body motions 
by the number of zero eigenvalues. In the 2-D case, there are 3 independent rigid body motions: 
translation in x- and y-axes and rotation. 

Discretization of the Displacement Field 


All the integral calculations at the element level (Fig.Al) are carried out in the isoparametric 
local coordinate system (f, 77 ), then transformed into the global cartesian coordinate system (x, y). 
For a function /, this transformation is: 

f fdV = t I f f(x,y)dxdy = t f f \ J \ d£dri . (16) 

Jv Jr a Jy a J-l J - 1 
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| 7 | is the determinant of the Jacobian matrix for isoparametric transformation, t is the thickness 
of the elastic plate. The problem is set up as C° and naturally considers only linear degrees of 
freedom (u,v). They are i-parallel and y-parallel displacements and may not be {-parallel and 77 - 
parallel. Plane linear isoparametric 4 -node quadrilateral serendipidity elements with bilinear shape 
functions are used to discretize the displacement field (Fig.A 2 ): 

Ni = j(l + &0(1 + ViV), (17) 

{,• = [- 111 - 1 ], 

77 , = [ -1 1 1 - 1 ], 
i = [ 1 2 3 4 ]. 

Any point in the element and its displacement can be related to the nodal point coordinates and 
their displacements such that: 


and alternatively, 
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(18) 


(19) 


( 20 ) 

( 21 ) 

( 22 ) 
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As noted in Eqn.8, the strain field is related to the displacement field by the differential operator 


D, 


D = 


a 

37 


a 

57 


0 

a 

37 

a 

37 


(23) 


The engineering definition for shear strain is implied; strain being: 
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(24) 


Because the calculations are carried out in the natural(isoparametric) coordinate system, the 
strains have to be transformed into this coordinate system. This can be established only by trans- 
forming derivatives of the displacement vector. However, the only way to transform derivatives of 
displacements is by applying the chain rule. Let 9 be a function of x and y, then invoking chain 
rule yields: 


(25) 


or 




y.z 


1 

H 

1 



x ,i y,v 


J 


u 

9„ 


= J 


And its the inverse relation is: 


',x 


= r 


u 

9.n 


(26) 


(27) 


where 


r = J~ l , 


(28) 


or 
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y,v -y,t 



(29) 


| J | is the determinant of the jacobian matrix, J. Depending on the nodal point coordinates 
assigned, the size, shape and orientation of the element will change, and so will the elements of the 
Jacobian. 

It is apparent that 0 is either u or v. The displacements in the natural coordinate system then 
are related as: 

Tu o o 

^21 T22 0 0 U,|J 

0 0 T2i Tn v t ( 

o o r 2 r 12 . . u„ 

After finding the displacement function in the natural coordinate system, the nodal point dis- 
placements are found to be: 

<7i 
<72 

u ,4 Ni ( 0 N 2( 0 N 3^ 0 N 4( 0 <73 

u , r, _ N itl 0 N 2 „ 0 N 3 „ 0 N 4r) 0 <7 4 

v, 4 0 N u < 0 0 N 3( 0 N 4( q 5 

. “v J L 0 N lt „ 0 N 2 „ 0 N 3 „ 0 N 4 „ J <7 6 

<71 

The formulation above, recalling Eqn.10, explicitly defines the displacement-strain matrix, B , 
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0 110 
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*2., 
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*3., 

0 
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( 32 ) 


Calculations of 7 £o , y„ 0 , and 7* are straightforward, however, x- and y-consistent forces from 7, 
have to be calculated at their nodes. The present study, does not include surface tractions. For con- 
tinuous media contact forces P are zero. (Note that P is another unknown as well as displacements). 


The integrals in the equations are calculated by Gaussian Quadrature (Dahlquist and Bjorck, 
1974). With this technique, a polynomial of degree (2n — 1) is integrated exactly by n-point 
Gaussian quadrature. The Legendre polynomials are used to solve for the coefficients( referred to 
as Gauss-Legendre coefficients). The sampling points and weight in Gauss-Legendre integration 
are tabulated(for example, see Bathe, 1982). For the element described above, a 2 by 2 point 
integration is sufficient. However, for very distorted elements, a 3 by 3 point integration might 
be preferred(Bathe, 1989; pers. comm.). Once the nodal point displacements are calculated, the 
strains for each element are obtained from : 


e = B q. (33) 

To calculate the stresses, one point integration for 4-node quadrilateral element at its center 
where £ = 7? = 0 is sufficient and most accurate! for arguments, see, for example Cook, 1981). The 
stress at the center of the element, from Eqn.l is : 


a = EBq-Ee 0 + cr 0 . (34) 

From the strain and stress tensors, principal major and minor strains and stresses, and their 
orientation can be calculated and will be used to complement the deformation pattern. Let \ denote 
the tensor whose principal components (when there is no shear) and orientation is desired. Using 
Mohr circle representation (Johnson, 1970), the tensor rotation yields the principal components 
*1,2> 


Ai,2 = ^(A,-A 2 )*±Ai, 


(35) 
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the maximum shear component A3, 


A 3 = \ II Aj — A 2 || , (36) 

and their orientation 0, 

0 = iarctan[j^-^-] . (37) 


The deformation pattern of any structure, hence, is displayed with the following plots: 

1. Undeformed structure with boundary conditions(scale in meters pointing North) 

2. Deformed structure, 

3. Displacement vectors(scale in meters) 

4. Principal major and minor) stress principal directions and their magnitude (scale in Pa), 

5. Magnitude of Maximum shear stress (scale in Pa), and 

6. Principal( major and minor) strain principal directions and their magnitude (scale has no 
units). 

Solution of Equations 

Follows a discussion of the solution of finite element equation(Eqn.lS) for P = 0, which corre- 
sponds to continuous media. For stable solutions, at least 3 degrees of freedom must be constrained. 
Also, for structural stability, prescription of certain other degrees of freedom, depending on the 
boundary/internal conditions might be required( directly affecting qr). The prescription of initial 
stresses (not included in the current study), body forces and initial strain also effect Q. Practically, 
the prescription of displacements at some nodes suppresses rigid body motion, as well as reducing 
the rank of the global stiffness matrix. At these nodes, the reaction forces must be calculated. 

For continuous media, Eqn.15 prevails. First, the solution techniques for continuum are dis- 
cussed, the deformation pattern for continuum exampled and, following the arguments favoring 
discontinuum approach , the solution technique and algorithms presented. 
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Consider Eqn.15 and partition the displacement vector q into two parts: q a and qt, representing 
the unknown and prescribed displacements, respectively. The nodal force vector is partitioned in 
the same sense: Q a and Qb, where the latter represents the unknown reaction forces. A similar 
partitioning of the stiffness matrix is also performed. This leads to: 


’ Kaa 

Kab ' 




’ Q*' 

m Kba 

Kbb _ 


. qb . 


_ Qb . 


The unknown displacements are 


(38) 


q a ~ R aa [Q a K a bqb] ? 

and the unknown reaction forces are 


(39) 


Qa = K a b q a + Kbb qb • (40) 

This is a direct solution for the unkowns and could be solved by any linear system solvers which 
utilize the bandwidth of the stiffness matrix. In this study, LU decomposition technique is used 
(Press et ah, 1986). However, when the structure gets complex, the the advantage provided by the 
bandwidth properties of the stiffness matrix is lost and makes direct solution techniques very costly. 
Therefore, an accelerated iterative solution becomes more attractive. Let <p define the residual of 
the / — th iteration: 


then, the ith degrees of freedom(displacement) at the Ith iteration is 


(41) 


= ?! _1 - « d , {Vi|i£Af}, (42) 

where n is the total number of degrees of freedom, A/* is the set of prescribed degrees of freedom, 
and a? is the acceleration factor. This is known as the successive over- relaxation method. Later, 
for the contact algorithm, <p will represent the unknown contact forces developing at the contacting 
nodes. In calculation of <p, the sparsity of the stiffness matrix is utilized to save computational 
time. Any mode in a finite element mesh connects a finite number of elements and as a result, the 
displacement at this node is affected by the nodes of these elements. So, the multiplication amount 
is reduced more than 60 % by skipping the nodes which are not connected. Theoretically, the 
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maximum number of iterations to converge and the related acceleration factor, can be calculated 
from the spectral radius of the relaxation operator. The large number of degrees of freedom yields 
(jj ~ 2.0, and this causes divergence in all solutions. The number of maximum iterations will also 
be underestimated. These parameters are problem dependent and u> ~ 1.6 (Cook, 19781). For the 
problems concerned, u> ~ 1.8 - 1.9. Cost is reduced at least 5- fold and the convergence which is 
controlled in each iteration by(for Eqn.15): 

c = £ II v‘i II • («) 

f 

C is the measure of convergence and must decrease as more iterations are performed. 
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Figure Captions 

Fig. 1: (a) The type of plate interactions at their boundaries in x-y-z space, (b) 2-D finite 
element approximation of plate boundaries when the dip of their boundary is 90° for divergent 
and trancurrent, and continent-continent collison boundaries , and 0 a for continent-ocean collison 
boundary [ x-z plane], (c) The behavior of plates at their boundaries in the x-y plane. Tension re- 
lease(seperation), sliding contact(slip), sticking con tact(single- node-continuum) and overlap modes 
corresponding to divergence, strike-slip faulting, continent-ocean collision and continent-continent 
collision, respectively. 

Fig. 2: The deformation pattern of a continuum model simulating transcurrent plate boundary. 

Fig. 3: The deformation pattern of a continuum model simulating divergent plate boundary. 

Fig. 4: The deformation pattern of a discontinuum model simulating transcurrent plate bound- 
ary: contact problem. 

Fig. 5: The deformation pattern of a discontinuum model simulating divergent plate boundary: 
contact problem. 

Fig. 6: Schematic representation of contact problem, after Bathe and Chaudry, (1985). 

Fig. 7: Coulomb Navier Criterion, after Wang and Voight(1969). 

Fig. 8: Nodal configuration at the contact surface 

Fig. 9: Deformation pattern of third type plate boundary: contact problem mode II. 

Fig.10 : Deformation pattern of second and third types plate boundaries: contact problem 
modes II and III. 

Fig. 11 : Deformation pattern of first type plate boundary: contact problem mode I. 

Fig. 12 : Deformation pattern of fourth type plate boundary: contact problem mode IV. 

Fig. 13 : Deformation pattern of fourth type plate boundary: contact problem mode II with 
internal deformation. 

Fig.14 : Deformation pattern of third type plate boundary: contact problem mode II with body 
forces caused by elevation changes. 

Fig.15 : Present-day tectonics of Eastern Mediterranean. After Hempton, 1985. 

Fig.16 : Focal Mechanism Solutions. After McKenzie(1978) 

Fig.17 : Local relative plate velocities. 

Fig. 18 : Index map of the region where the finite element contact problem is applied. 

Fig. 19 : Deformation pattern of the eastern Mediterranean[Model 1]: Contact. Driving force: 
boundary displacements. 
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Fig.20 : Deformation pattern of the eastern Mediterranean[Model 2]: Contact. Driving force: 
boundary displacements. The convergence at the Cyprean arc and shortening at the Bitlis Suture 
are included. 

Fig.21 : Deformation pattern of the eastern Mediterranean[Model 3]: Contact. Driving forces: 
boundary displacements and gravitational forces. The convergence at the Cyprean arc and short- 
ening at the Bitlis Suture, internal deformation at the Palmyra Kink are included. 

Fig.Al: Global and isoparametric(natural) coordinate systems and their mapping. 

Fig.A2: Bilinear shape functions for 4-node quadrilateral element. 
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Fig. 1: (a) The type of plate interactions at their boundaries in x-y-z space, (b) 2-D finite 
element approximation of plate boundaries when the dip of their boundary is 90° for divergent 


and trancurrent, and continent-continent collison boundaries , and 0° for continent-ocean collison 
boundary [ x-z plane], (c) The behavior of plates at their boundaries in the x-y plane. Tension re- 
lease( seperation), sliding contact(slip), sticking contact(single- node-continuum) and overlap modes 
corresponding to divergence, strike-slip faulting, continent-ocean collision and continent-continent 


collision, respectively. 
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Fig. 4: The deformation pattern of a discontinuum model simulating transcurrent plate bound- 
ary: contact problem. 
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Fig. 5: The deformation pattern of a discontinuum model simulating divergent plate boundary: 
contact problem. 
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(c) Forces acting on contactor and target bodies 
Figure 1. Schematic representation of problem considered 


Fig. 6: Schematic representation of contact problem, after Bathe and Chaudry, (1985). 




after Wang and Voight(1969). 




Fig. 8: Nodal configuration at the contact surface 
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Fig. 9: Deformation pattern of third type plate boundary: contact problem mode II. 
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tig.10 : Deformation pattern of second and third types plate boundaries: 
modes II and III. 
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Fig. 11 . Deformation pattern of first type plate boundary: contact problem mode I. 
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fig. 12 : Deformation pattern of fourth type plate boundary: contact problem mode IV 
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Fig. 14 : Deformation pattern of third type plate boundary: contact problem mode II with body 
forces caused by elevation changes. 
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Fig- 15 : Present-day tectonics of Eastern Mediterranean. After Hempton, 1985. 
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Fig. 17 : Local relative plate velocities. 



tig. 18 : Index map of the region where the finite element contact problem is applied. 
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tig. 19 : Deformation pattern of the eastern Mediterranean[Model 1]: Contact. Driving force: 
boundary displacements. 
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Fig.20 : Deformation pattern of the eastern Mediterranean! Model 2]: Contact. Driving force: 
boundary displacements. The convergence at the Cyprean arc and shortening at the Bitlis Suture 
are included. 


ORIGINAL PAGE IS 
OF POOR QUALITY 



"** ' - 1 • *0>. rt Uj«l*«ti ,..i.,| .... .« , 4 


»ti Sh««. Sir. 


• - >■* t*jin 



tig. 21 : Deformation pattern of the eastern Mediterranean(Model 3]: Contact. Driving forces: 
boundary displacements and gravitational forces. The convergence at the Cyprean arc and short- 
ening at the Bitlis Suture, internal deformation at the Palmyra Kink are included. 
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Zone 

A. Aykut Barka and M. Nafi Toksoz 

Earth Resources Laboratory 
Department of Earth and Planetary Sciences 
Massachusetts Institute of Technology 
Cambridge, MA 02142, USA. 


Abstract 

Historical and instrumental earthquakes in the eastern part of the North Anatolian 
fault zone between the Erzincan basin and the Karliova triple junction have been 
examined in relation to fault segmentation and kinematics. The 12/26/1939 Erzincan 
earthquake (M=8) created a 360 km surface break and it was terminated at the 
eastern end of the Erzincan pull-apart basin in the east. The 8/17/1949 (M=6.< - 
7) earthquake was a double bend earthquake which affected the easternmost three 
fault segments (FS1, FS2 and FS3) of the North Anatolian fault zone. According 
to the historical data, the 1784 earthquake (I=VIII-IX) occurred between the 1939 
and 1949 rupture segments (FS4-FS9, about 75 km long). From this data it appears 
that the North Anatolian fault zone could have two seperate sequences of westward 
migrating large earthquakes. One of the two sequences originates from the tripie 
junction between the NAFZ and the Ovacik fault in eastern end of the Erzincan basin, 
and extends to the west about 900 km, as happened between 1939-1967. The second 
one extends between the Karliova triple junction and the Erzincan basin and consist' 
of two rupture segments, 1949 and 1784. The main reasons for these two separate 
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migrations are firstly different recurrence intervals of large earthquakes of the rupture 
segments which is controlled by the geometry and length of the rupture segment and 
secondly, the eastern part of the westward escaping Anatolian block divided into two 
wedge shaped blocks (A1 and A2) each of which moves independently. From historical 
records it may appear that the recurrence interval of the western migration could be 
about 900 years while the recurence interval is about 200-300 years for the eastern 
migration. Furthermore, the eastern migration has not yet completed during the 
current twenteeth century migration. Thus, within the eastern sequence, the 1784 
rupture segment to the east of the Erzincan basin is identified as a potential seismic 
gap. Recurrence intervals of historical earthquakes and geological data indicate that 
slip rate in this part of the fault zone is about 0.8-1 cm/yr. This results in an 
accumulation of about 2 m right-lateral slip along the 1784 rupture segment. 

Introduction 

It has recently been emphasized that fault geometry plays a critical role in the earth- 
quake rupture process (e.g., Segall and Pollard, 1980, Bakun et al. 1980, Lindh and 
Boore 1981, Barka and Hancock 1982, King and Nabelek 1985, Sibson 1986, Schwartz 
and Coppersmith 1986, Barka and Kadinsky-Cade 1988, Wesnousky 1989). The term 
’’fault geometry” includes stepovers, bends, and their many combinations. In this 
study we focus on strike-slip fault geometry and earthquake activity in the eastern 
part of the North Anatolian fault zone. We have studied the geometry of active fault 
segments in detail on the field and aerial pothographs, belonging not only the North 
Anatolian fault also other major faults in the region. Then, we combined this infer 
mation with distribution of intrumental and historical earthquakes in the region. We 
also examined the extents of the surface ruptures of the large earthquakes thorough 
the litrature and some on the field. As a result our purpose is to understand how the 
fault moves and examine seismic gaps along the fault zones. 

Figure 1 shows major tectonic elements of Turkey in an area where the northward 
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motion of the Arabian plate causes active convergence in eastern Turkey. As a result, 
the Anatolian block escapes westward and the Northeast Anatolian block eastward 
(Ketin 1948, McKenzie 1972, Sengor 1979, Kasapoglu k Toksoz 1983, Gulen 1984, 
Dewey et al. 1986). The wedge shaped Anatolian block is bounded by the right- 
lateral North Anatolian fault to the north, and by its conjugate, the East Anatolian 
fault, to the south. These two fault zones intersect at the Karliova Triple junction 
(Ketin 1969, Allen 1969, 1975, McKenzie 1972, Dewey 1976, Tchalenko 1977, Sengor 
1979, Toksoz et al. 1979, Jackson k McKenzie 1984, Sengor et al. 1985, Dewey et al. 
1986). The eastern part of the Anatolian block is divided into two smaller blocks ( Al 
and A2 in Fig. 2) by the left-lateral strike-slip Ovacik fault. This fault intersects the 
NAF zone at the southeast end of the Erzincan basin which forms an an other triple 
junction (ETJ1). The eastward escape of the NE Anatolian block is complicated by 
the extensive internal deformation and by the existence of a number of sub-blocks. 
The Northeast Anatolian fault zone (NEAFZ), forms the northern boundary of the 
NE Anatolian block. The dominant tectonic feature in this region is the NAF, which is 
a joint boundary between the two blocks escaping in opposite directions. The NAFZ 
intersects the NEAFZ in northwest of Erzincan (ETJ2, Figures 1 and 2). Figure 
2 shows the major blocks and boundary faults in the area of concern, between the 
Erzincan and Karliova triple junctions. Genuinly, both historically and during the 
modern times, the Erzincan region has been one of the most active seismic regions in 
Turkey (Sieberg 1932, Ergin et al. 1967, Soysal et al. 1981, and see Table 1) because 
the area is situated within a most critical tectonic center from where continental 
blocks escape sideways. 

Between 1939 and 1967 most of the North Anatolian Fault west of Erzincan rup- 
tured through a westward migrating series of major earthquakes, as shown in Figure 
1. In this series of earthquakes, the largest one was the 1939 Erzincan earthquake 
(M=8.0). East of Erzincan, earthquakes along the NAF followed a more complicated 
pattern, as can be seen in Figure 2. 
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Fault and Rupture Segments 


Fault segmentation is defined by the distribution of the geometric discontinuities along 
straight fault segments. For the definition of minumum sizes of these discontinuities 
which control the fault segmentation, we used Barka and Kadinsky-Cade (1988) crite- 
ria (stepover widths and bend angles larger than 1 km and 5 respectively). Rupture 
segments are the extents of surface rupture zones produced by characteristic large 
earthquakes. The North Anatolian fault zone consists of a number of fault and rup- 
ture segments in this area, as shown in Figure 2 (Barka and Kadinsky-Cade, 1988). 
In this section we outline the geologic and seismic details of each fault and rupture 
segment belonging to the major fault zones. 

The North Anatolian Fault Zone 

The NAFZ initiates at the Karliova triple junction (Figs. 1 and 2). Although there is 
extensive seismic activity to the east of this junction (between Karliova and Varto). 
it is believed that this part is no longer the continuity of the NAFZ, but a rather 
complex suture zone which has developed by the westward escape of the Anatolian 
block. In other words, one interprate this a is the triple junction was initially in the 
Varto area and as a result of the westward escape the Anatolian block, the triple 
junction has moved to the Karliova area. Barka and Gulen (1988) have called this 
zone as a “zipper zone”. The complex suture zone has been formed by a thrust 
formation extending through the bisector of the angle between the boundary strike- 
slip faults (Fig. 3) which is the manifastation of closure of the space at the tip of 
the escaping block. The 1946 and 1966 Varto earthquakes occurred along this zone 
(Tasman 1946, Ambraseys and Zatopek 1968, Wallace 1968, Ketin 1969 and see Figs 
2 and 3). According to McKenzie (1972) and Canitez (1973) the fault plane solution of 
the 1966 Varto earthquake differed from other solutions of eartquakes which occurred 
along the NAFZ by having a trust component with right-lateral strike-slip motion 
(Fig. 4). This is in good agreement with suture zone formation. In the same context. 
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an E-W trending thrust morphology was interprate from aerial photos in the Karliova 
area extending through the bisector angle between the boundary strike-slip faults, 
NAFZ and EAFZ, (Fig. 3). This north dipping thrust seems to splay from the 
segment 1 of NAFZ and indicates that the Karliova basin can not be an extensional 
basin, instead it is a complex ramp basin formed by the interaction of coeval strike- 
slip faulting and thrusting. Although the epicenter of 7/07/195 < earthquake, M=5.5, 
is located 30-50 km west of the Karliova area, its thrust mechanism (Canitez and 
Ucer 1967), indicates that the internal deformation of the Anatolian block includes 
an approximately E-W thrusting which can be used as supporting evidence for the 
thrust formation in the Karliova area (Fig. 4). 

Fault segments 1, 2 and 3 (FSl, FS2 and FS3) form a restraining doublebend 
between the Yedisu basin and Karliova triple junction. Double bend angles are 20 
in east and 25 in the west (Fig. 3). FSl extends from the triple junction to the 
west about 30 km and has very clear physiographic expressions. Along this segment 
many small streams and ridges which are normal to the fault trace are offset and 
curved in a right-lateral sense (see also Allen 1969 and 1975). This segment, the 
fault trace is expressed by a long norrow trough along its entire lenght (Fig. 3l 
FS2, which is the restraining segment of the double bend, has two small releasing 
stepovers and runs within the Elmali river. Its physiographic expressions are much 
less clearer developed. FS3 forms the western part of the double bend and also has a 
very clear strike-slip morphology like FSl (Fig. 3). The 8/20/1966 Ms=5. 3-6.2 and 
some other smaller aftershocks of the 1966 Varto earthquake were located on h! 
(Ambrasesy and Zatopek 1968, Dewey 1976). During the 8/20/1966 aftershock mo>t 
of the villages in the vicinity of this segment were destroyed (Fig. 3). The fault plan*- 
solution of this earthquake (Mckenzie 1972, Canitez 1973) indicated a pure strike -dtp 
motion along FSl. The 8/17/1949 Ms=6.7-7 earthquake affected all three segment 
According to data collected during our field survey and Ambraseys (1987. person.*! 
communications) this earthquake might have created ruptures mostly along FS2. ami 
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the eastern part of FS3 and the western part of FSl. Ambraseys (1987, personal 
communation) also reported that right-lateral displacements might have reached up 
to 1.5-2 m in places. Most of the damage occurred along FS2 which makes up the 
restraining area of the doublebend (Fig. 3). The fault plane solution of this event 
had a slight thrust component (Canitez 1973). 

FS3 and FS4 form the 2.5 km wide Yedisu pull-apart basin (Fig. 5). There is a 
restraining angle about 10 between these two segments. FS4 is about 28 km long and 
has clear morphological expresions. FS4 has also two small stepovers first of which is 
a releasing type located in the Yedisu basin and the second one is a restraining type 
situated in middle of the segment. The 7/26/1967 Pulumur earthquake, M=5.6-6.2, 
took place west of this second stepover and created 4 km surface breaks and 20 cm 
right-lateral displacement along this part of FS4 (Ambraseys 1975). Figure 5 shows 
intensity contours and destroyed villages along the FS4 (Tutuncu and Derrurtasli 
1967). The fault plane solution of this earthquake (Mckenzie 1972 and Canitez 19/3) 
indicated that the motion was also pure strike-slip and confirms the right- lateral 
motion of the NAFZ (Fig. 4). 

FS4, FS5 and FS6 form a releasing double bend (Fig. 5). In the vicinity of 
FS6 there are several other faults which trend parallel to FS6. There is a long and 
narrow small lake along FS5 which is consistent with its extensional nature. FS7. 
FS8 and FS9 create a combination of releasing bend and restraining stepover. The 
width of the restraining stepover is about 2 km. FS7, FS8 and FS9 have relatively 
less clear morpological expressions. This is probably due to fact that they all run 
within the Ephratus valley where fluvial activity is quite high. Fault expressions 
are well developed only between Tanyeri and Caykomu villages along FS9 (Fig 5i 
Ambraseys (1975) reported that the 1784 large earthquake occurred to the east of the 
Erzincan basin and created 90 km long surface rupture. Based on this information it 
is believed that this rupture segment extended along the entire lenght of FS4-F c ''t 
In other words, 1784 rupture segment took place between the Yedisu and Erzinran 
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releasing stepovers. 

FS9 and FS10 form a releasing bend (15) and stepover (4-5 km) combination 
which is responsible for the opening of the Erzincan basin (Fig. 2). FS10 is 60 km 
long, and has clear physiographic expressions in its western half. The southeast half 
of FS10 is characterized by short en-echelon strike-slip faults and contemporaneous 
volcanics (Barka and Gulen 1989). FS10 is separated from FS11 by a 20 restraining 
bend. FS11 is about 110 km long, and extends from this bend, situated about 10 
km NW of the Erzincan basin, to Susehri - the location of a releasing double bend 
(Kocyigit 1988). The area of interest in this paper terminates in the eastern part of 
FS11 (for the details of the other fault segments to the west see Barka and Kadinsky- 
Cade 1988). FS11 has clear strike-slip morphology especialy in the vicinity of Mihar 
village. Along all these fault segments the strike-slip motion is associated with a 
vertical component. According to field observations, the southern block is usually 
uplifted except that this varies where the fault segments form extensional structures 
(releasing bend and/or releasing stepover). 

The 1939 Erzincan earthquake created surface ruptures along FS10, FSll and 
extended further west along the Kelkit valley and terminated south of Amasya where 
there is a 24 restraining bend between the fault segments (e.g. Pamir and Ketin 1941. 
Parejas et al. 1941, Ketin 1969 see also Figure 1). The total lenght of surface rup- 
tures was about 360 km. Although, Parejas et al. (1941) reported 3.7 m maximum 
displacement, according to recent studies (Kocyigit 1988, Barka in prep.) the max- 
imum displacement reached 7.5 m in along FSll. During the same earthquake the 
southern block was uplifted about 2 m . Pamir and Ketin (1941) reported that two 
foreshocks were felt within the week preceding the main shock in the Erzincan basin 
The epicenter of the earthquake was located near the 20 restraining bend, on FS10 
(Dewey 1976). Many of the 1939 earthquake aftershocks caused damage in the Er.i- 
incan and Niksar pull-apart basins (Nature 1940a, b, c, d, Ergin et al., 1967; Tabban. 
1980; see also Riad and Meyer, 1983). Fault plane solution of this earthquake 
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pure strike-slip motion and agreed with the motion along the fault zone. An another 
fault plane solution for a moderate size earthquake (Mb = 4.8, 11/18/1983) near the 
city of Erzincan is characterized by ENE-WSW extension (International Seismologi- 
cal Centre Bulletin solution 1983), also in agreement with the pull-apart character of 
the Erzincan basin (Fig. 4). 

The East Anatolian fault zone 

The left-lateral East Anatolian fault zone is the southern boundary of the westward 
escaping Anatolian A1 block. The most northeastern fault segment (FSl) is about 
60 km long and extends from the Karliova triple junction to Bingol (Figs, land 2). 
The segment is straight and has clear morphological left-lateral strike-slip expressions 
and it is also accompanied by a normal component (western block down) along the 
northern half. The southern half of this segment runs within the Goynuk river valley 
where the expressions are not so clear. The 1971 Bingol earthquake (M=6.7) created 
surface ruptures mostly along the southern half of the segment (Arpat and Sarogiu 
1972 and Seymen and Aydin 1972). The fault plane solution of this earthquake 
indicated pure left-lateral slip along this segment (Fig. 4). At least one other historical 
earthquake (1789) of a similar size has been documented from Soysal et al. ( 198 1 ) in 
the vicinity of the same segment. According to Arpat and Sarogiu (1972) and Seymen 
and Sarogiu (1972) the fault zone has 15-27 km left-lateral post-Miocene displacement 
revealing about 0.5 cm/yr slip rate. 

The Ovacik fault 

This is another left-lateral fault and is about 120 km long trending NE-SW. According 
to Barka and Gulen (1989) who studied the tectonic evolution of the Erzincan basin, 
the Ovacik fault has also been participating in the opening of the Erzincan basin 
The Ovacik fault splays into several branches before it enters the Erzincan basin 
The Ovacik basin is situated on a releasing bend along the segments of the Ovacik 
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fault. In the Ovacik basin the fault cuts Quaternary alluvial fans, and forms very 
distinct 5-10 m high fault scarps (see also Arpat and Saroglu, 1975). Although the 
fault, in general, has a left-lateral strike-slip character, these scarps indicate that the 
motion in the Ovacik basin is also associated with a normal component. However, 
outside the Ovacik basin the fault has a thrust component that causes the uplift of 
the Munzur Montains. As far as historical earthquakes are concerned there are no 
spesific events in the last 1000 years that can be associated with this fault. During the 
present field survey occurence of one large earthquake 1200 years ago was interpreted 
from the Legend of the Munzur Springs in the Ovacik area. Barka and Gulen (19S9) 
estimated 5-7 km left-lateral displacement along the Ovacik fault from the geomet ry 
of the southeastern part of the Erzincan basin. They also considered that the age of 
the Ovacik fault is younger than the NAFZ (approximately 3-3.5 Ma). These values 
may reveal about 0.15-0.25 cm/yr slip rate along the Ovacik fault. 

The Northeast Anatolian Fault Zone 

This fault zone consists of several segments with a combined length of approximately 
350 km. The southwesternmost segment (FS A) is located to the north of the Erzincan 
region (Figure 2). Approximately 70 km long, it strikes NE-SW. Although very little is 
known about this fault segment, it is assumed to have an oblique movement, consisting 
mostly of left-lateral slip with a subordinate thrust component (Tatar, 1978). The 
study of earthquake records (Soysal et al., 1981; Sipahioglu, 1983; Riad and Meyers. 
1985) indicates that it is less active than the segments of the North Anatolian Fault 
zone. Pamir and Ketin (1941) showed ESE-WNW trending isoseismals parallel to 
the NAFZ, covering the area between Tercan and Baskoy for the 11/21/1939 Tercan 
earthquake (M=5.9). And because of this, this earthquake was always considered 
to be the foreshock of the 1939 Erzincan earthquake. However, after locating the 
damaged villages and the main trace of the NEAFZ, we believed that this earthquake 
may have occurred on FS-A of the NEAFZ and has no relation with the NAFZ < l*»r 
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example, 130 buildings collapsed in Karakulak which is situated next to the fault 
zone; Pamir and Ketin, 1941; Ergin et ah, 1967; Tabban, 1980, see also Figure 2 
where the locations of other destroyed villages are shown). This is also confirmed by 
the relocation of the earthquake (Dewey 1976, Fig. 7). Apart from the 1939 Tercan 
earthquake and several aftershocks of the 1939 great Erzincan earthquake, the only 
known historical event associated with this segment is the 1254 earthquake ( 1= IX). 
This event caused surface breaks to occur over a 50 km length on FS-A (Ambrasevs, 
1975). 


Seismicity 

Historical Earthquake Records 

The history of damaging earthquakes in the Erzincan region was recognized and 
fairly well documented even before the great earthquake of 1939 (Ali Kemal, 1932). 
Sieberg (1932) listed some of the Erzincan earthquakes and stated that between 104 5 
and 1784, at least 17 damaging earthquakes had occurred in the Erzincan region In 
Table 1 we have tabulated the significant earthquakes affecting the Erzincan region 
since 1000 A.D., based on sources referenced in the table. Figure 6 is an intensity-time 
plot of known earthquakes which have affected the Erzincan region. From this figure, 
earthquakes can be categorized according to two large sizes: (a) great earthquakes 
for which I X (Modified Mercalli intensity), and (b) large earthquakes with VIII 1 IX 
According to Figure 6, at least 3 great earthquakes have occurred during the last 10u0 
years, including the one in 1939. Ambraseys (1970) reported that the 1045 earthquake 
produced a surface break of length comparable to the one which occurred in 1939: and 
that the 1458 earthquake caused the death of about 32,000 people, however, although 
this figure is comparable to the casualties of the 1939 earthquake, in most catol»<g> 
the affected area is described as taking place between Erzincan and Erzurum. 

At least 10 large earthquakes (VIII I IX) have occurred in the Erzincan rcg i.:. 
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since 1000 A.D., causing considerable damage and large numbers of casualties. Among 
the earthquakes of this size only 1254 and 1784 earthquakes can be associated with 
specific segments; 1254 was on FS-A of the NEAFZ and 1784 was on FS4-FS9 of 
the NAFZ (Ambraseys 1975). The other two large earthquakes similar 1784 were 
1578 and 1422 which were separeted by about 156 years. Furthermore, there is the 
possibility of only one large earthquake along the Ovacik fault occuring about 1200 
years ago which is interpreted from the Legend of Munzur Springs. In other words, 
none of the listed earthquakes in Table 1 is specifically associated with the Ovacik 
fault. 

The recurrence interval (900 year) combined with displacement created during 
the great earthquakes (7.5 m), give a slip-rate of approximately 0.8 cm/yr. This slip- 
rate is similar to that obtained from geological data which reveals about 0.8—1 cm/\r 
for this part of the fault zone (Barka and Gulen 1988). Note also that FS1-FS10 
form a joint boundary between opposite-moving blocks (the Anatolian and Northeast 
Anatolian blocks). Thus the slip rate is expected to be higher in this area than along 
the main section of the NAF to the west. However, some amount of the slip should be 
taken up by the internal deformation of A1 block as expressed by extensive internal 
faulting and related seismic activity at tip of the wedge shape block. From Figure 
6 the recurrence interval for large earthquakes (VIII I IX) is approximately 100-150 
years. 


Instrumental Earthquake Records 

Figure 7 shows the distribution of earthquakes M 4.9 in the region between 1900 and 
1985. Epicenters for the interval 1900-1930 are taken from Tabban (1980) and Riad 
and Meyers (1985). Epicenters of those earthquakes which occurred between 1930 
and 1971 are taken from Dewey (1976) who relocated the events M 5. Moreover, 
epicenters of all earthquakes for the period of 1964-1984 belonging to ISC, are also 
shown in Figure 8. From both maps two significant points can be made; a) marn 
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moderate- large earthquakes are mostly concentrated along the NAFZ, and b) there is 
also concentrations of the small- moderate earthquakes at the tips of the wedge shaped 
blocks. The tip of A1 block has the most clear activity out of the escaping blocks. 

Seismic Gaps 

Both historical data and the 1939 earthquake have shown that great earthquakes in 
this region appear to be consistent with the 1939 rupture segment which includes 
FS10, FSll and fault segments (FS11-FS14) to the west along the NAFZ (see Barka 
and Kadinsky-Cade 1988). If we consider the recurrence interval of 1939 earthquake 
(M=8) to be about 900 years (Fig. 6) this rupture segment is safe for long time for 
a similar size of earthquake. On the other hand there is no historical data obtained 
before the 1949 earthquake which took place at the eastern end of the fault zone 
(Fig. 9). This is propably due to that the area is sparsly populated because of its 
raged morphology. However, if we take the estimated slip rate as about 1 cm/yr and 
similar amount of slip (1.5-2 m) on the fault, one can simply calculate 150-200 years 
recurence interval for the 1949 size of earthquakes. Along the 1 (84 rupture segment, 
before 1784 earthquake, there are two other comparable earthquakes which affected 
the Erzincan region within last 500 years( 1578, I = VIII and 1422, 1= VI II, Fig. 9). but 
their locations are not known. The recurrence interval of these earthquakes are 156 
and 206 years. Since it has been 205 years since 1784 earthquake, this rupture segment 
(FS4-FS9), stands out as a clear seismic gap. The geometry of the fault segments also 
indicate that the restraining features along the rupture segment are not large enough 
to restore so much stresses. Thus, one can conclude that the geometry and historical 
earthquake records and long term slip rate along the rupture segment indicate that 
this gap should have a large earthquake near future. The estimated 1 cm/yr slip rate 
results in over 2 m slip accumulation along this rupture segment. An other important 
point with this gap is that during the 20th century this rupture segment is the only 
segment along the NAF zone which has not experienced a large earthquake between 
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Varto and the western end of the Mudurnu valley (900 km). Moreover, this gap along 
the NAFZ is different from the gap mentioned by Toksoz et al., (1979), and was first 
briefly mentioned by Ambraseys and Zatopek (1969). 

There is no any significant earthquake that can be specifically associated with the 
Ovacik fault since 1900 during the instrumental period and/or last 1000 years. Thus, 
it has been 1200 years since the last known large earthquake along this fault. With 
the given rate (0.15-6.25 cm/yr), about 1. 5-2.5 m left slip might have accumulated 
along this fault. Therefore, the Ovacik fault may well be another candidate for future 
large earthquakes. 

The amount of total displacement along the FS-A of the NEAF zone is similar to 
the Ovacik fault (about 5 km). The 11/21/1939 Tercan earthquake and 02/03/1949 
aftershock of the 1939 Erzincan earthquake might have occurred on this segment. 
From the historical earthquake records, we are only aware of the 1254 large earth- 
quake, which created 50 km of surface faulting along segment A, trending 60 with 
5 m (?) maximum vertical displacement (Ambraseys, 1975). The slip rate with the 
historical data indicates that about 1 m left-lateral slip may have been accumulated 
along this rupture segment. 

Migration of Large Earthquakes 

The historical data is not long enough to clearly understand the migration patterns of 
the North Anatolian fault zone. However, from available data, two separate westward 
migration of large earthquakes along the North Anatolian fault can be integrated 
One starts from the triple junction of the North Anatolian and the Ovacik fault 
(ETJ1) and to the westward as it happened between 1939 and 1967. The second 
sequence takes place between the Karliova and Erzincan triple juction (ETJl) which 
consists of two rupture segments 1949 and 1784. The western migration might occur 
approximately every 900 years while the eastern migration may repeat every 200 100 
years. The longer recurrence interval of 1939 rupture segment is related to the 20 


13 



restraining bend which is the largest restraining feature along the North Anatolian. 
Furthermore, two separate Anatolian blocks (A1,A2) move to west thus two separate 
sets of migration earthquakes would be expected when each block moves. As far as 
most recent sequences concerned as has mentioned the eastern migration has not yet 
completed. 

Fault Geometry and Earthquake Rupture Processes 

Some of the details of the fault geometry and rupture processes have been already 
discussed by Barka and Hancock 1982, Barka and Kadinsky-Cade 1988 and Kadinsk\ - 
Cade and Barka 1989. Some of these can be sammurized as follows, 

a) It is apparent that two ends of the 1784 rupture segment is controled by the 
Erzincan and Yedisu pull-aparts. 

b) Each rupture segment has restraining feature in itself, the size of which propor- 
tional to the size of earthquake, length of rupture segment and the amount of slip For 
example, as has mentioned the above, the characteristic great Erzincan earthquakes 
(1045 and 1939) are closely related to the 20 restraining bend (110 km long) north- 
west of the Erzincan basin and the restraining bend of the 1939 is larger than those 
observed along the 1949 (20 and 15 km long) and 1784 rupture segments so as the 
size of earthquakes, the rupture lengths and the displacements. More details of these 
issues have been discussed by Barka and Kadinsky-Cade (1988), Kadinsky-Cade and 
Barka (1988) and Wesnousky (1988). 

c) As it has been also poited out that location of epicenters, in other words, rupt ure 
initiation take place mostly nearby the locked segments (e.g. Barka and Hancock 
1982, King and Nabalek 1985, Barka and Kadinsky-Cade 1988). For example, the 
epicenter of the 1939 earthquake is located near the restraining bend (Barka and 
Hancock, 1982; Barka and Kadinsky-Cade, 1988). Similarly, the epicenter of I'M'* 
double bend earthquake is located on FS3, near the the western bend where angle 
higher relative to the eastern bend. Thus, from these examples one can assume ihu’ 
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for the case of the 1784 seismic gap the epicenter of the large earthquakes may take 
place in western part of FS4 and/or FS5 (Fig. 9). 

d) Furthermore, as mentioned earlier, in the easternmost part of the fault zone 
FS1, FS3 and FS4 have clear morphological expressions. Many stream beds and 
ridges normal to the segments are curved in a right-lateral sense. These three seg- 
ments have a linear geometry and also are parallel to the slip direction. It should be 
noted that most of the area of interest including the Ovacik fault and segment A of the 
NEAF zone, is located within the serpentinite-rich ophiolites and ophiolitic melange 
associated with the Anatolide/Tauride - Pontide suture zone. These three aspects, 
morphological expressions, linearity in geometry and plasticity of the rocks indicate 
high potential of creep phenomena along these segments. In other words the motion 
along the straight segments is considered to be easy. The seismicity indicate that 
this assumed creep is accompanied with continuous small and moderate earthquake 
activity, if we consider that those moderate earthquakes (such as the 1966, M=6.3 
along the FS1 and the 1967 Pulumur earthquake, M=5.6-6.2, along the FS4, Figure 
9) are characteristic earthquakes of these straight segments . Note that for any case 
these activities (creep and small-moderate earthquakes) do not exclude the potential 
for future large earthquakes along the FS1-FS9 as it was experienced in the 1949 and 
1784 events. However, we believe that the amount of slip during those large earth- 
quakes should be less along the those segments that have creep and small- moderate 
earthquake activity. 

In summary, one can postulate that straight segments which are parallel to the 
slip direction has interseismic activity of small- moderate earthquakes and creep, and 
they continuously transfer the slip (or stresses) to the locked areas (Fig. 9). Locked 
segments who has restraining geometries move coseismicly by large events ( F i g ‘ ■ 
During the large earthquake those straight segments which have already transfer the 
slip to the locked areas would have obviously less slip. For example, if the above 
mention gap along the North Anatolian fault creates a large earthquake, we would 
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expect maximum slip to the west of FS4 and much less slip along the FS4. 

Finally, we can conclude that releasing stepovers or releasing double bends appear 
to be site of preseismic and post seismic activity such as foreshocks and aftershocks 
(Fig. 9). If this is so, the 1784 seismic gap could have a high potential of foreshock 
activity and most of the aftershock could also take place in the Erzincan and Yedisu 
pull-aparts. 


Conclusion 

From this study it is clear that examination of earthquakes and fault geometry pro- 
vide many useful information in understanding not only the tectonics of the side way 
escapes of the continental blocks but also to in defining seismic gaps along the major 
fault zones which form the boundaries of blocks and earthquake rupture processes. 

This study also indicates that there is a clear 75 km long seismic gap along the 
North Anatolian fault immidiately to the east of the Erzincan basin. This segment 
last ruptured in 1784 and since the estimated slip rate is about 1 cm/yr, about 2 m 
slip has accumulated. The recurence interval of large earthquakes along this rupture 
segment is an average of 180 years. Although the damage and casualties were less 
severe than 1939, the 1784 earthquake was also very destructive for the Erzincan 
region killing at least 5 000 people. The other significant seismic gap appears to be 
associated with the Ovacik fault along which 1. 5-2.5 m slip might have accumulated 
over last 1200 years. 

There could be two separate westward migration of large earthquakes along the 
North Anatolian fault zone, first one is from the Erzincan basin to the west and second 
one is from Karliova triple junction to the Erzincan basin. This is caused by primarily 
different recurrence intervals of large earthquakes along the rupture segments and two 
separate block motion to the west. 

Slips along the straight segments which are also parallel to the general slip di- 
rection seem to be easy and expressed by clear morphological expressions, frequi-:.' 
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small- moderate earthquakes and possible creep activity trnsfering the stresses to the 
locked segments. Finally, a high potential of foreshock activity is associated with the 
pull-apart structures and a releasing double bend along the the gap segment. 
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Figure Captions 


Figure 1. Tectonic map of Turkey showing the surface ruptures due to major earth- 
quakes since 1900. The Anatolian and Northeast Anatolian Blocks are wedged out 
to the west and east respectively by the convergense of Arabia and Eurasia as shown 
in the inset map (lower left). The rectanle in the Figure delinates the area of study 
and is enlarged in Figure 2 (compiled from Arpat and Saroglu 1972, 1975, Arpat et 
al. 1977, Barka and Hancock 1984, Sengor et al. 1985). 

Figure 2. Simplified geometry of major blocks and distribution of fault and rupture 
segments between Erzincan and Karliova. Thick and dashed lines and dates indicate 
rupture segments and dates of related earthquakes, respectively. Doted areas are 
Plio-Quaternary basins being formed mostly by strike-slip faulting. Al and A2 are 
sub-blocks within the Anatolian block. Stars near the rs 1939b show the destroyed 
villages by the 11/21/1939 earthquake (taken from Ergin et al. 1967). 

Figure 3. Simplified tectonic map of the Karliova- Varto area, rs 1949, fs2 are 
rupture segments of major earthquakes and fault segments, respectively (see also 
Tutkun 1986 and Saroglu et al. 1987). The observed extents of surface ruptures of 
the 08/19/1966 Varto earthquake is indicated by dashed lines (from Wallace 19681 
The triangles are the villages destroyed by the largest aftershock of the 1966 Varto 
earthquake. Closed squares within the squares are the villages knocked down by the 
08/17/1949 earthquake (information collected during the present survey and Am 
braseys 1987 written communication). 

Figure 4. Fault plane solution of major earthquakes of the region (compiled from 
Mckenzie 1972 and Canitez 1973). 

Figure 5. Characteristic features of the North Anatolian fault between Yedisu and 
Erzincan basins. Doted contours are isoseismals of the 07/26/1967 Pulumur earth 
quake and stars are the destroyed villages both taken from Tutuncu and Derrurta.-h 
(1967). Squares are the destroyed villages by the 08/17/1949 earthquake. 

Figure 6. Time (T)/intensity (I) distribution of earthquakes in the Erzincan area 
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Numbers above the dotes are the number of casualties resulting from each particular 
event, a and b are the categories of earthquakes. For explanation and references see 
the text and Table 1. 

Figure 7. Distribution of earthquake epicenters (M > 4.9) in the easternmost part 
of the North Anatolian fault zone for the interval 1900-1987. Solid circles indicate 
the epicenters taken from Dewey (1976) and the open circles indicate epicenters taken 
from mostly Riad and Meyers (1985) and Tabban (1980). 

Figure 8. Distribution of ISC epicenters of all earthquakes between 1964-1984 in 
the region. 

Figure 9. Summary of the relationship between geometry of the fault segments and 
seismicity along the eastern part of the North Anatolian fault zone. Horizontal axis 
represents the extent of the fault and vertical axes are time and slip. The continuous 
vertical lines indicates each rupture segments and dashed vertical lines separate the 
straight segments and locked segments. Stars with dates ilustrate locations relocated 
epicenter of the 1939 and 1949 earthquakes (Dewey 1976). The star with exclamation 
mark is the location of epicenter of the expected large earthquake within the gap area 
The horizontal arrows indicate the direction of the slip by moderate earthquakes. The 
dates, size and the amount of slip of the large earthquakes are also indicated. 
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Table 1. List of historical earthquakes 
in the Erzmcan Region. 


Number Date Intensity (I) Number of casualties 


(1) 

1045 

X-Xl 


(2) 

1161 

VI 


(3) 

1165 

VII 


(4) 

1166 

VI 

12.000 

(5) 

1160 

VIII 

(8) 

1170 

viii-ix 


(7) 

1236 

VI 


(8) 

1251 

VIII 


(9) 

1254-55 

VIII 

16,000 

(10) 

1268 

IX 

15.000 

(U) 

1207 

VIII 


(12) 

1209 

VIII 


(13) 

1308 

VI 


(14) 

1358 

V 


(15) 

1366 

VI 


(18) 

1374 

VII 


(17) 

1422 

VIII 


(18) 

1433 

VI 


(19) 

1458 

X 

32,000 

(20) 

1543 

VII 


(21) 

1579 

VIII 

1,500-15.000 

(22) 

1605 



(23) 

1667-9 

vm-x 

Half of the town 



was destroyed 

(24) 

l 7 94 

VIII-IX 

5.000-15.000 

(25) 

1997 

VI 



* Docxr.e.T.ed ?to.t. : 332 V. '.932 So oxoi-Civ. 1936-1940, Pa_-«jas «£ al., 1941 

Pinar and Lahn 1952, Zrgsi «i at ,96^ W:-a3eyi 1970 1975, <attlx 1972. Can 1974. 

Soysad^Xoi., 1961 ! 962 5 . : an. -g . ^ '967 
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Slip Distribution of the Great 1939 Erzincan 
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The December 26, 1939 Erzincan earthquake (M=8) is the first and largest of a 
remarkable westward migrating series of six M=7-8 earthquakes that occurred along 
the North Anatolian fault zone between 1939 and 1967 (Ketin, 1948, 1969; Ambraseys 
1970). The 1939 earthquake created a 360 km long surface break, which is shown in 
Figure 1 (Ketin, 1969; Barka et al., 1987). The epicenter was located approximately 
10 km northwest of Erzincan, as indicated by the star in Figure lc (Dewey, 1976). The 
distribution of aftershocks is not known completely, although a number of moderate 
size aftershocks can be associated with the Erzincan, Su§ehri and Niksar basins, 
which are major releasing features along the rupture zone (Dewey, 1976; Tabban, 
1980; Riad and Meyers, 1985). Several small earthquakes were felt by local residents 
in the Erzincan basin two weeks before the main shock (Pamir and Ketin, 1941). Only 
one measurement of fault slip was made soon after the earthquake: a 3.i m right- 
lateral offset of a side wall along the main road in Re§adiye (Parejas et ah, 1942), 
In this short note we present results from recent field measurements of fault offset 
along the 1939 earthquake surface breaks. These data are very important, not only 
because they help constrain the rupture characteristics of the 1939 earthquake, but 
also because they allow us to compare the 1939 earthquake with other great strike-slip 
events such as the 1857 and 1906 San Andreas fault earthquakes in California. 
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The rupture zone can be divided into five major fault segments (Barka and 
Kadinsky-Cade, 1988). This division is illustrated in Figure lc. From east to west 
the lengths and strikes (measured from north) of these segments are. (1) Erzincan 
length 60 km, strike 125°, (2) Mihar - length 65 km, strike 105 , (3) Ortakoy - length 
45 km, strike 117°, (4) Kelkit - length 100 km, strike 107°, and (5) Niksar - length 90 
km. strike 90°. These five segments are related geometrically in the following fashion: 
(a) segments 1 and 2 are separated by a 20° restraining bend (Barka and Hancock, 
1982), (b) segments 2, 3 and 4 form a releasing double bend along which 3 is the re- 
leasing segment, and (c) segments 4 and 5 are separated by a 17° smooth restraining 
bend. 

During the field survey we measured offsets of man-made features such as field 
boundaries, fences, roads, canals and lines of trees defining field boundaries. We also 
measured offsets of natural features such as streams and the side walls of valleys 
or ridges. The survey focused on villages that were situated right on the surface 
breaks. Thus local residents who had witnessed the earthquake were able to confirm 
our identification of surface breaks as well as some of the offset features. 

The 1939 earthquake slip distribution determined from the new measurements is 
summarized in Figure lb. Both man-made and geomorphological offset features are 
included in the figure. An independent set of measurements of surface slip associated 
with the 1939 earthquake is provided by Kocyigit (1988). These data points are 
included in Figure lb as well. They are in close agreement with our measurements. 
From the combined data set shown in Figure lb, the surface slip distribution can be 
described in the following simplified fashion: 4 - 6 m slip along the western half of 
segment T, increasing to 6.5 - 7.5 m along segments 1 and 2, then decreasing to 4 - 
4.5 m along segment 4, followed by a further decrease to 2 - 2.5 m along segment 5 
(possibly less at the western end of segment 5). 

Figure la shows a reversed seismogram from the 1939 earthquake (time shown 
increasing from right to left for purposes of discussion). The seismogram was recorded 
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in Pasadena, California. The epicentral distance for this record is 104°, so that 
the phase shown in the figure is a diffracted P wave. The semi-major axis of the 
90% confidence ellipse for Dewey’s (1976) relocated epicenter coordinates (39.80 N, 
39.38°E) was only 10-20 km, so we are quite certain that the epicenter can be tied to 
t lie 20° restraining bend just northwest of Erzincan (Figure lc). The Pasadena record 
starts with an approximately 20 second long low-amplitude wavetrain (emergent phase 
starting at the ”1” label in Figure la and lasting until ”2”). Most of the moment 
release, however, occurs during the next 100 seconds (starting shortly after 2 and 
lasting until ”3”), with the largest amplitudes occurring in the first 60 seconds of that 
100 second wavetrain. The simplest way to interpret this seismogram is to associate 
the 20 second long low-ampltude phase with rupture of the 60 km segment (1) located 
near Erzincan, and the 100 second long phase with rupture of the remainder of the 
fault zone (west of the epicenter). This interpretation is based on (1) a reasonable 
rupture velocity of 3 km/sec (60 km in 20 seconds for the low-amplitude phase and 
300 km in 100 seconds for the higher-amplitude phase), and (2) a good correlation 
between the distribution of amplitudes in the seismogram and the distribution of 
surface slip (Figures la and lb). Clearly, however, a more thorough study of historical 
seismograms produced by this earthquake needs to be done. 

The slip distribution along the fault segments suggests that maximum slip is 
associated with the restraining (west) side of the 20° restraining bend. That side 
of the bend is uplifted, and folding and thrusting are common features in the late 
Cenozoic sediments of that area (Tatar, 1978; Barka and Gulen, 1989). The relocated 
epicenter near the bend supports a model of bilateral rupture propagation. The fault 
plane solution of the earthquake determined by McKenzie (1972) is characterized 
by predominantly right-lateral strike-slip motion with a small component of reverse 
faulting. The strike of the main fault plane in his solution is 108°, similar to the strike 
of segment 2 in Figure lc. 

The distribution of slip shown in Figure lb can be utilised to estimate a static 
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moment for the 1939 earthquake, using the formula Mo=/iuA, where /x is the rigidity 
of the medium, u is the average dislocation and A is the area of faulting. We assume a 
rigidity of 3.3 x 10 n dyne/cm 2 and a crustal thickness of 15 km, and add the moments 
from each of the 5 fault segments. Here average slips of 5 m, 7 m, 7 m, 4.25 m and 
2.25 m are assumed for segments 1 through 5 based on Figure lb because a more 
complicated calculation is not warranted by the data. 

Mo = (3.3 x 10 n dyne/cm 2 )(15Jfcm)[(60km x 5m) + ((65 + 45)km x 7m) + (100km x 
4.25m) + (90km x 2.25m)] = 8.4 x 10 27 dyne.cm 

The corresponding moment magnitude based on the formula 

log Mo = l.bMw +16.1 

of Hanks and Kanamori (1979) is Mw=7.9. 

The slip deficit to the west along segment 4 and 5 can be interpreted in one or 
both of the following ways, (a) It is possible that another significant earthquake has 
occurred along this section of the fault in the past, or could occur in the future, to 
make up the deficit, (b) Internal deformation of the Anatolian Block may be taking 
place, and the slip deficit may correspond to deformation occurring at the nearby 
Tokat kink. As supporting evidence for the latter posibility Barka and Gulen (1988) 
have pointed out that total displacement along the fault zone decreases from 35-40 km 
in the Erzincan region to 25-30 km in the central section of the North Anatolian fault. 
However it should be noted that the region of slip deficit along segment 4 coincides 
with the eastern portion of the surface break produced during the 1668 earthquake 
(Ambraseys and Finkel 1988). 

There is a fundamental difference between the slip distribution of the 1939 earth- 
quake and those of the 1857 and 1906 earthquakes along the San Andreas fault. In 
the case of the California earthquakes it is possible to compare the fault geometry 
with slip distributions reported by Thatcher (1975) and Sieh (1978). In the case of 
the 1857 earthquake the section of the fault between Hwy. 166 and Tejon Pass can 
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be described as the restraining section of a double restraining bend. Here slip was 
only about 6 m, compared with 9 m in the Carrizo plain to the north. In the case 
of the 1906 earthquake a similar situation occurred in the southern section of the 
rupture zone, just north of San Juan Bautista. Reduced slip in the restraining sec- 
tions suggests that these areas acted as barriers to rupture propagation. In the 1939 
F.rzincan earthquake case the restraining section just west of the epicenter was the 
area of maximum slip. That section of the fault may have been acting as an asperity 
in 1939 and as a barrier in 1668. 
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Figure Caption 


Figure 1. Rupture characteristics of the 1939 Erzincan earthquake, (a) Reversed 
(time increasing right to left) version of Pasadena (PAS) record of diffracted P from 
the 1939 earthquake (Z component, To= 1 sec, Tg=0.23 sec, time correction= -20 
sec). The origin time of the earthquake was December 26, 23:57:16.0 (GMT). For 
explanation of numbers above seismogram see text, (b) Slip distribution of the 1939 
M = 8 Erzincan earthquake. Solid circles correspond to slip measured during this 
survey. Open squares are measurements of Kocyigit (1988). The open triangle is 
the only measurement obtained soon after the earthquake (Parejas et ah, 1942). (c) 
Major fault segments of the 1939 rupture. The inset map shows the location of the 
1939 rupture zone relative to the trace of the North Anatolian fault. 
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APPENDIX 4 


PRELIMINARY REPORT ON 1989 TURKEY GPS CAMPAIGN 

This is a very brief, preliminary report on the results of the MIT 1989 Turkey GPS 
field campaign. In short, due to the exceptional effort of UNAVCO personnel, the excel- 
lent cooperation and logistical support provided by the Turkish Union of Geodesy and 
Geophysics (TUJJB), and the dedicated effort of the individual field parties, the mea- 
surement campaign was extremely successful. In spite of the failure of one of the Trimble 
receivers, we were able to observe all 18 stations planned for this year’s campaign. This 
was accomplished by extending the survey by 4 days - a possibility which arose when the 
Greek/ Aegean experiment was postponed. 

The accompanying Table lists sites observed this year by MIT/TUJJB. Continuous ob- 
servations were made at ANKARA (ANKA) and DIYARBAKIR (DIYA) using TI-4100 
receivers, while the other sites in Eastern Turkey were observed for 3 days each using 
Trimble 4000ST receivers. A total of 4 station-days of data were lost, out of 86, due to 
instrument and logistical problems. As these losses occurred on 3 separate days and on 
4 different stations, they will not significantly reduce the strength of the network. 

Our GPS measurements were closely coordinated with those made by IFAG and their 
collaborators at and around SLR sites in Turkey, as well as with the measurements in 
Western Turkey made by the Durham University group. The accompanying Figure shows 
the locations of stations observed by our group in 1988 and 1989, those observed by IFAG 
and 5 of the approximately 30 sites observed by Durham (these 5 sites were established 
and observed with GPS by MIT/Hacettepe Univ. in 1988). As indicated in the Figure, 
the 1989 sites established by MIT/TUJJB in Eastern Turkey are well located to monitor 
regional deformation associated with on-going continental collision in this area. This 
includes: 

1. The distribution of crustal shortening between the Arabian plate and the Eurasian 
plate along a transect running from the Turkish/Syrian border to the Black Sea, 

2. “Extrusion” and rotation of the Anatolian plate with concentrated deformation 
along the North and East Anatolian faults, 

3. “Extrusion” rotation and internal deformation of the East Anatolian block, and 

4. Relative movement between the African and Arabian plates along the Dead Sea 
fault. 

In addition, strong ties were established to our 1988 GPS observations in Western Turkey 
through overlapping observations at 5 sites in central and Eastern Turkey (4 SLR sites 
and Ankara) and the 5 1988 GPS stations observed this year by Durham. Furthermore, 
continuous observations at the SLR stations in Askites (ASKI) and Rhodes (RllOD) 


made by IFAG will provide ties to the Aegean network. 

Besides completing a very successful measurement campaign, we established a close work- 
ing relationship with the Turkish geodetic community. The 1989 campaign was done in 
close cooperation with the TUJJB. This is an important development for the project 
both for assuring the long term viability of our measurement program in Turkey and for 
receiving increased logistical support for future surveys. We feel strongly that the success 
of this project, as well as all international projects, requires more of the field effort be 
taken on by local scientists. The relationships developed this year are an important step 
in this direction. 

We will begin processing and analyzing the Turkey GPS measurements using the GAMIT 
software this month. One Turkish scientist, to be selected by TUJJB and MIT, will visit 
MIT to participate in data reduction. One of our first efforts will be to compare GPS 
and SLR baselines, as all 4 SLR sites in Turkey were reobserved by both techniques this 
year. We are also very interested in comparing 1988 and 1989 GPS observations at SLR 
sites, the Ankara site, and the 5 GPS sites in Western Turkey observed in 1988 and 1989. 
We expect that these studies will be undertaken in cooperation with the other groups 
involved with GPS and SLR measurements in this region. 
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Figure 1. 


Some of the GPS sites observed in and around Turkey in 
Key: • 1989 MTT/TUJJB; Q 1988 MIT/Hacettene U./IFAG; 


1988 and 1989..-. 
• 1989 IF AG ; (J 


1989 Durham U. 


SLR stations in Rhodes (RHOD) and Askites (ASKI) were observed by J^G in 1988. 
nth#>r 1988 sites were observed bv MIT/Hacettepe . See Table for details o 
SS/™” sUel? Snly 5 of approximately 30 Dufham sites are shown (these 5 srtes 
were established and observed by MIT/Hacettepe^n 1988> IFAG sites clos ^ 

stations (Footprints) shown schematicallv. DION, ASKI, RHOD, hua, 

DIYA are SLR sites. 


